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Abstract 

Optical  Tomography  of  the  ZaP 
Flow  Z-Pinch  Plasma 

by  Adam  Matthew  Madson 

Chair  of  the  Supervisory  Committee: 

Associate  Professor  Uri  Shumlak 
Department  of  Aeronautics  &  Astronautics 

The  ZaP  Flow  Z-Pinch  Experiment  forms  a  Z-pinch  with  an  embedded  axial  flow 
that  exhibits  stability  for  500-2000  times  the  instabililty  growth  time  for  a  static 
Z-pinch.  Two  linear  photodiode  arrays  placed  orthogonal  to  each  other  measure 
chord-integrated  emissivity  of  the  plasma.  Deconvolving  these  data  produce  time 
resolved  1-D  and  2-D  emission  profiles  of  the  Z-pinch  cross-section  during  the  sta¬ 
ble  quiescent  period.  These  unfiltered  emission  profiles  show  that  the  pinch  has  a 
radius  of  approximately  1  cm.  Filtering  the  light  for  Ha  emission  shows  a  defined 
hollow  profile  1  cm  in  radius.  Bremsstrahlung  filtered  light  exhibits  a  peaked  profile 
during  periods  of  high  intensity  emissivity  and  a  hollow  profile  at  lower  intensities. 
Additional  measurements  from  an  Imacon  camera,  interferometer,  and  azimuthal  ar¬ 
rays  of  magnetic  probes  determine  the  spatial  correlation  between  plasma  emissivity, 
plasma  density,  and  plasma  current.  The  tomographic  information  collected  by  the 
photodiode  arrays  increases  understanding  of  flow  stabilization  theory  and  Z-pinch 
configuration  plasmas. 
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Chapter  1 

INTRODUCTION 

With  the  decline  of  traditional  fossil  fuel  reserves  around  the  globe,  there  is  great 
interest  in  the  research  and  development  of  new  “alternative”  energy.  A  safe,  clean, 
reliable  energy  with  a  nearly  unlimited  or  even  renewable  source  is  the  objective  of 
current  research  in  the  energetics  held.  Fusion  reactions,  the  processes  where  the 
combination  of  light  atomic  nuclei  into  heavier  nuclei  release  substantial  amounts 
of  energy,  offer  an  exciting  area  of  research  in  this  held.  These  reactions  naturally 
occur  inside  the  sun,  where  predominately  hydrogen  ions  form  helium  nuclei  through 
the  proton  chain  process  [11].  Unfortunately,  this  is  a  slow  and  weakly  interacting 
process,  and  it  is  only  the  massive  central  pressure  inside  the  sun  that  increases  the 
reaction  rate  to  a  level  that  sustains  fusion.  The  gravity  created  by  the  sun  contains 
this  reaction,  providing  a  counter-balance  to  the  outward  force  resulting  from  the 
pressure  gradient  between  the  ionized  gas  at  the  core  and  at  the  surface.  Therefore, 
a  successful  fusion  reaction  depends  on  the  containment  of  a  highly  dense  volume 
of  mostly  ionized  gas,  called  plasma.  As  defined  by  Francis  Chen,  “A  plasma  is  a 
quasineutral  gas  of  charged  and  neutral  particles  which  exhibits  collective  behavior 
[3].”  Several  different  containment  schemes  for  highly  dense  plasmas  are  currently 
under  investigation  around  the  world,  including  Spheromak,  Tokamak,  Reversed  Field 
Pinch,  and  Z-pinch  plasmas.  The  ZaP  experiment  at  the  University  of  Washington 
Aerospace  and  Energetics  Lab  attempts  to  contain  high  density  Z-Pinch  plasmas 
using  an  innovative  sheared  axial  flow  configuration.  Two  multi-channel  photodiode 
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arrays  measure  optical  emissions  from  the  plasma,  comparing  them  to  other  diagnostic 
systems  as  well  as  providing  one-dimensional  and  two-dimensional  deconvolved  radial 
profiles. 


1.1  Introduction  to  the  Z-Pinch 


A  Z-pinch  is  a  relatively  simple  plasma  configuration  in  terms  of  shape  and  con¬ 
tainment  requirements.  It  is  a  cylindrically  symmetric  plasma  column  with  an  axial 
current,  Jz,  as  shown  in  Figure  1.1.  This  axial  current  induces  an  azimuthal  mag¬ 
netic  field  Bg  which  in  turn  confines  the  plasma  radially.  Without  any  plasma  and 
resulting  current,  there  would  be  no  magnetic  field.  Therefore,  the  Z-pinch  is  entirely 
self-contained,  with  no  externally  driven  magnetic  fields  needed  for  its  creation.  The 
simplicity  of  this  configuration  greatly  reduces  the  complexity  of  the  force  balance 
equations  governing  the  equilibrium  pressure  balance. 

Ampere’s  law  and  the  ideal  magneto-hydrodynamic  (MHD)  force  balance  equation 
are 

VxB  =  fi0J  (1.1) 

J  x  B  =  Vp.  (1.2) 


Given  that  Jz  and  Bg  are  the  only  non-zero  current  and  magnetic  field,  these  equations 
reduce  respectively  to  the  following  [6] 


Jz  =  -4-  (r Bg) 
Hqt  dr 

T  p  dP 
dr 


Substituting  1.3  into  1.4  yields 


d_  (  B$_\  _  B$_ 

dr  V  2/h)/  IMyr ' 


(1.3) 

(1.4) 


(1.5) 


The  derivative  terms  represent  the  plasma  particle  pressure  and  the  magnetic  pressure 
respectively.  The  plasma  pressure  force  is  always  pointing  out  from  the  radial  axis 
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Figure  1.1:  Z-pinch  geometry 


of  the  Z-pinch.  The  next  term  is  the  magnetic  pressure  force,  which  points  outwards 
when  confining  the  plasma  outside  the  plasma  radius.  The  term  on  the  right  hand 
side  is  the  magnetic  tension  term,  created  by  the  curvature  of  the  magnetic  field 
lines.  For  a  discrete  plasma  with  a  uniform  current  profile  out  to  its  radius,  the 
plasma  and  magnetic  pressure  forces  balances  the  magnetic  tension  force  outside  the 
plasma  column.  An  easy  way  to  illustrate  magnetic  tension  term  is  by  visualizing  the 
force  generated  by  a  rubber  band  wrapped  around  an  object.  It  is  this  inward  force 
that  provides  the  radial  force  balance  needed  to  confine  a  Z-pinch.  This  balance  also 
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results  in  greater  confinement  efficiency,  represented  by  the  term  (3. 

0  average  plasma  pressure  <  p  > 
peak  magnetic  pressure  Bq/ 2p0 

This  confinement  makes  the  Z-pinch  a  good  candidate  for  fusion  power  applications, 
unfortunately,  real  life  instabilities  also  make  it  violently  unstable. 


1.2  Instabilities  Inherent  to  a  Z-pinch 

Z-pinches  are  unstable  to  both  current  driven  kinks  and  pressure  driven  MHD  insta¬ 
bilities  [6].  These  instabilities  arise  when  a  displacement  in  the  plasma  offers  a  lower 
energy  state  versus  the  existing  configuration.  Any  variations  in  the  current  profile  or 
the  Bq  held  can  lead  to  instability.  The  two  best  examples  are  the  m  —  1  and  m  —  0 
instability  modes. 

The  m  =  0  mode,  or  sausage  mode  shown  in  Figure  1.2(a),  occurs  when  a  small 
section  of  the  plasma  column  radially  contracts.  The  reduced  cross-sectional  area 
increases  the  magnetic  field  at  this  location  and  therefore  its  magnetic  tension  force 
directed  inward.  This  results  in  even  further  narrowing  of  the  plasma  at  this  location, 
until  the  plasma  physically  breaks  apart. 

The  m  —  1  mode,  or  kink  mode  shown  in  Figure  1.2(b),  occurs  when  the  plasma 
column  bends,  creating  a  stronger  local  magnetic  held  on  the  inner  side  of  the  bend 
and  a  weaker  held  on  the  outer  side.  The  resulting  increase  in  magnetic  tension 
creates  a  force  in  the  direction  of  the  kink,  further  increasing  its  deformation  un¬ 
til  containment  is  no  longer  possible  [6].  Traditional  methods  for  stabilization  have 
included  modifying  the  pressure  profile,  reducing  wall  radius,  and  applying  axial  mag¬ 
netic  fields.  Unfortunately,  these  methods  all  have  significant  drawbacks  that  reduce 
their  usefulness  for  fusion  plasma  applications.  Limiting  the  radial  pressure  gradient 
for  example,  can  stabilize  the  plasma  against  the  rn  —  0  mode.  Specifically,  if  the 


(a) 


(b) 


Figure  1.2:  (a)  rn  —  0  “sausage”  mode  and  (b)  rn  —  1  “kink”  mode  with  surrounding 
Be  field 


pressure  profile  satisfies  the  following  relationship 


4r  d\np 

2  +  T(3  ^  d  In  r 


(1.7) 


where  T  is  the  ratio  of  specific  heats  and  j3  is  the  ratio  of  plasma  pressure  to  magnetic 
pressure,  the  m  =  0  mode  can  be  marginally  stabilized  [15,  19].  Unfortunately,  this 
has  no  effect  on  stabilizing  the  m  =  1  mode. 

Likewise,  an  imbedded  axial  magnetic  field  in  the  plasma  can  counteract  the  mag¬ 
netic  tension  of  the  m  =  1  mode,  but  it  also  limits  the  maximum  axial  current  and 
opens  up  the  field  lines  [16,  17].  As  the  plasma  column  bends,  the  axial  field  stretches, 
providing  a  larger  and  larger  force  to  counteract  the  instability.  Limiting  of  the  cur¬ 
rent  lowers  the  maximum  attainable  plasma  pressure,  making  this  option  undesirable 
for  fusion  power  applications.  These  drawbacks  in  conventional  stabilization  methods 
resulted  in  the  reduction  of  Z-pinch  research  in  favor  of  other  plasma  configurations. 
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1.3  Sheared  Flow  Stabilization 

Theoretical  research  by  Shumlak  and  Hartman  [18]  and  experimental  research  by 
Shumlak  et  al.  [19],  provides  evidence  that  sheared  velocity  flow  offers  the  exciting 
possibility  of  limiting  instability  propagation  in  the  Z-pinch.  Specifically,  a  radial 
shear  of  the  axial  flow,  dvz/dr  >  O.l/cVd,  where  a  is  the  plasma  radius,  k  the  wave 
vector,  and  Va  the  Alfven  velocity,  can  satisfy  the  marginally  stable  pressure  profile 
condition  previously  given  by  Equation  1.7.  Alfven  velocity  relates  to  the  electro¬ 
magnetic  Alfven  wave.  An  out  of  plane  electric  field  produced  by  the  wave  combines 
with  the  Bo  magnetic  field  to  produce  a n  E  x  B  drift  [3] .  This  results  in  the  plasma 
particles  moving  in  the  same  oscillating  variation  as  the  magnetic  field  lines,  motion 
known  as  frozen-in  flow.  Perturbation  propagation  occurs  as  these  variations  in  the 
plasma  travel  down  the  field  lines  and  cause  local  instability  modes.  The  speed  of  this 
propagation  is  given  by  the  Alfven  velocity,  defined  by  Va  =  B0/ y/pop  [3].  There¬ 
fore,  relating  the  Alfven  speed  to  radial  velocity  shear  allows  for  the  construction  of  a 
normalized  stability  criterion  for  varying  radii.  This  criterion  is  shown  in  Figure  1.3. 


1.4  Overview  of  the  ZaP  Flow  Z-pinch 

The  ZaP  Flow  Z-Pinch  Experiment  at  the  University  of  Washington  Aerospace  and 
Energetics  Research  Laboratory  currently  investigates  the  effects  of  sheared  axial 
flow  on  Z-pinch  plasma  stability.  In  its  present  design,  a  coaxial  acceleration  region 
similar  to  a  Marshall  gun  in  an  attached  plasma  assembly  region  produces  1  meter 
long  Z-pinch  plasma  columns.  The  present  design  of  the  experiment  is  included  in 
Figure  1.4 

To  form  the  Z-pinch,  hydrogen  gas  is  injected  into  the  vacuum  chamber  between 
the  inner  and  outer  electrodes  near  the  midplane  of  the  acceleration  region.  Applying 
a  negative  voltage  to  the  inner  electrode  causes  the  gas  to  ionize  and  form  a  current 
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Threshold  Shear 


Figure  1.3:  Normalized  velocity  shear  threshold.  rw/a  is  wall  radius  over  plasma 
radius,  v'  is  dv~/dr,  Va  is  the  Alfven  velocity,  and  k  =  2tt/\  is  the  wave  number. 
Plasmas  whose  edges  are  close  to  the  wall  are  generally  stable,  however  for  fusion 
applications  the  pinch  radius  is  much  smaller  than  the  wall  radius.  The  threshold 
stability  for  this  case  would  correspond  to  the  right  side  of  the  figure,  with  a  value  of 
0.1. 
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Figure  1.4:  The  ZaP  Flow  Z-pinch  Experiment 


sheet  running  from  the  outer  to  the  inner  electrode.  This  current  induces  a  magnetic 
held  in  the  azimuthal  direction,  and  the  two  combine  to  produce  a  J  x  B  Lorentz 
force  in  the  z  direction  that  pushes  the  plasma  down  the  inner  electrode  and  into 
the  assembly  region.  As  the  current  moves  into  this  region,  the  increasing  distance 
and  slant  between  the  inner  electrode  nosecone  and  outer  electrode  connection  point 
produces  a  component  of  the  Lorentz  force  pointing  towards  the  2- axis,  causing  the 
current  sheet  to  collapse  down  onto  axis.  Once  there,  the  plasma  pressure  increases 
and  the  Z-pinch  forms  with  an  imbedded  sheared  axial  how.  This  sequence  is  shown 
in  Figure  1.5. 

Neutral  gas  remaining  in  the  acceleration  region  continues  to  ionize  and  feed  into 
the  pinch  until  no  more  remains  or  instabilities  destroy  the  pinch.  Experimental 
results  show  that  the  ZaP  sheared  how  Z-pinch  remains  stable  500-2000  times  longer 
than  the  instability  growth  time  of  a  static  Z-pinch.  Other  characteristics  of  the  ZaP 
Z-pinch  are  included  in  the  Table  1.1  [14]. 
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Figure  1.5:  Z-pinch  formation,  (a)  Neutral  gas  is  injected  in  the  acceleration  region 
and  (b)  ionized  when  a  potential  is  applied  between  the  electrodes,  (c)  The  J  x  B 
Lorentz  force  accelerates  the  current  sheet  into  the  assembly  region  (d  and  e)  and 
as  it  extends  in  length,  the  resulting  inward  Lorentz  force  pushes  the  plasma  onto 
axis,  (f)  This  creates  a  dense  Z-pinch  column  aligned  on  axis  with  an  embedded  axial 
shear. 
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Table  1.1:  ZaP  Flow  Z-pinch  typical  operating  parameters. 


parameter 

symbol 

value 

ion  velocity 

^ion 

5-10  cm/ /is 

Alfven  speed 

VA 

12-15  cm/ /is 

quiescent  period  length 

T? 

20-40  /is 

total  temperature 

Te  +  Ti 

150-200  eV 

electron  number  density 

Tie 

1016-1017  cm“3 

edge  magnetic  field 

-®edge 

1.0-1. 8  T 

Z-pinch  radius 

a 

1  cm 

Z-pinch  length 

L 

1  m 

peak  current 

-^p  max 

150-250  kA 

peak  power 

S p  max 

0.2-1. 5  GW 

capacitor  charge  voltage 

Vmb 

5-9  kV 
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Chapter  2 

DIAGNOSTIC  SYSTEMS  AND  COMPARISONS 

The  ZaP  experiment  currently  uses  a  variety  of  diagnostic  systems  to  investigate  Z- 
pinch  plasma  stability.  These  include  measurements  of  the  magnetic  field  and  plasma 
current  using  axial  and  azimuthal  arrays  of  magnetic  probes,  plasma  currents  using 
Rogowski  probes,  chord-integrated  plasma  density  from  a  two  chord  interferometer, 
and  light  intensity  or  emissivity  from  CCD  and  ICCD  spectrometers,  photomultiplier 
tube,  fast  frame  camera,  and  photodiode  arrays.  Much  time  and  research  goes  into 
the  secondary  confirmation  of  data  taken  by  a  single  diagnostic  device.  Relating 
photodiode  emissivity  data  to  the  current  location  measured  by  the  magnetic  probes, 
the  electron  density  measured  by  the  interferometer,  and  the  emission  profile  taken 
by  the  fast  frame  camera  provides  insight  to  the  stabilization  theory  of  the  sheared 
flow  Z-pinch. 

2.1  Photodiode  Array  (PDA) 

Photodiodes  are  relatively  simple  devices  that  prove  valuable  in  converting  light  inten¬ 
sity  to  a  current.  A  photodiode  crystal  consists  of  a  surface  P-layer  and  a  subsurface 
N-layer.  Photons  striking  the  crystal  stimulate  electrons  in  these  layers  out  of  their 
energy  levels.  The  electron-hole  pairs  move  in  opposite  directions,  with  the  electrons 
toward  the  N-layer  and  holes  toward  the  positive  P-layer.  When  placed  in  a  closed 
circuit,  the  voltage  difference  caused  by  this  movement  creates  a  current  that  varies 
with  light  intensity.  As  the  resistance  of  this  circuit  approaches  zero,  the  current 
varies  more  and  more  linearly  with  the  changing  light  level,  as  shown  in  Figure  2.1 
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Negative  Bias 


Figure  2.1:  Photodiode  current  vs.  voltage  curve.  In  a  closed  circuit  with  no  resis¬ 
tance,  current  changes  linearly  with  respect  to  light  intensity. 


[10] .  The  depletion  layer  created  on  the  anode  side  of  the  photodiode  between  the  P- 
and  N-layers  slow  the  response  time  between  light  intensity  variations  and  the  cur¬ 
rent.  Introducing  a  reverse  bias  current  on  the  photodiode  replenishes  this  depletion 
layer,  thereby  reducing  signal  response  time.  More  discussion  regarding  the  biasing 
circuit  is  in  Section  4.1. 


2.2  Magnetic  Probe  Arrays  and  m  —  1  Mode  Measurements 


The  magnetic  probes  in  the  ZaP  experiment  measure  mode  activity  and  current  prop¬ 
erties  of  the  Z-pinch  by  measuring  the  changes  in  the  local  magnetic  field.  Using  the 
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relationship  described  by  Faraday’s  law  [9] 

I  E-  c/1  =  -  J  B  da,  (2.1) 

coils  of  wire  imbedded  in  the  outer  electrode  generate  a  voltage  based  upon  the  varying 
magnetic  field.  A  single  axial  array  of  probes  measures  the  current  sheet  distribution 
and  location.  Four  azimuthal  arrays,  three  in  the  assembly  region  with  eight  equally 
spaced  probes  and  one  in  the  acceleration  region  with  only  six  probes  measure  the 
azimuthal  mode  activity  [8].  Normalizing  the  m  —  1  magnetic  mode  by  the  aver¬ 
age  magnetic  field  at  the  axial  location  of  the  array  indicates  the  displacement  of 
the  Z-pinch.  This  is  especially  important  since  it  allows  for  the  creation  of  an  em¬ 
pirical  boundary  condition  between  a  stable  pinch,  defined  as  the  quiescent  period, 
and  unstable  regions,  defined  as  the  formation  and  disruption  phases.  The  MAT- 
LAB  code  magneticmode  .m  used  to  determine  magnetic  mode  activity  is  included  in 
Appendix  A. 3. 

2.3  Interferometry  and  Chord-Integrated  Electron  Density 

The  two  chord  Mach-Zehnder  configuration  interferometer  determines  electron  density 
by  measuring  the  index  of  refraction  with  two  coherent  laser  light  sources.  One  beam 
travels  through  the  plasma  and  then  recombines  with  a  reference  beam,  resulting 
in  constructive  and  destructive  interference.  Photodiodes  record  these  changes,  and 
since  the  phase  of  the  laser  light  shifts  proportionally  to  the  chord-integrated  electron 
density,  multiple  beams  traveling  through  the  plasma  at  a  single  axial  location  allow 
for  the  determination  of  the  line-integrated  density  through  different  radial  locations 
of  the  plasma. 

2.4  Imacon  Camera 

Taking  optical  measurements  of  the  plasma  is  simply  a  matter  of  recording  the  image 
viewed  through  a  window  in  the  vacuum  chamber  onto  some  sort  of  medium.  High 
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Figure  2.2:  Normalized  m  —  1  mode  displacement.  Quiescence  occurs  when  the 
displacement  is  below  the  empirical  value  of  0.2  as  shown  by  the  dashed  black  and 
green  lines.  This  corresponds  to  a  radial  displacement  in  the  tank  of  approximately 
1  cm  off  the  z-axis.  In  this  example,  Z-pinch  quiescence  occurs  from  37.0  —  85.4  /is. 


speed  photography  offers  an  excellent  means  of  accomplishing  this  task  and  storing 
a  great  deal  of  visual  data  on  a  small  section  of  him.  An  Imacon  Model  790  camera 
takes  images  and  stores  them  on  Polaroid  667,  3000  speed  him.  Data  are  stored  as 
a  series  of  framed  images  at  stepped  time  increments  or  as  a  streak  photo  along  the 
length  of  the  Polaroid  him.  An  electrical  circuit  in  the  Imacon  controls  the  geometric 
placement  of  the  image  on  the  him  through  the  use  of  “shutter”  and  “compensating” 
capacitor  plates.  For  a  framing  picture,  the  camera  electronically  places  images  at 
separate  locations  on  the  him.  In  the  case  of  the  streak  photo,  a  narrow  slit  is  used  to 
view  only  a  thin  vertical  section  of  the  plasma.  The  camera  then  sweeps  this  image 
along  the  length  of  the  him  at  a  known  linear  rate.  This  creates  an  image  of  plasma 
movement  in  the  vertical  direction  versus  time  in  the  horizontal  direction  of  the  him. 
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Chapter  3 

TOMOGRAPHIC  THEORY 

The  use  of  photodiodes  to  collect  radiated  data  is  just  one  of  a  variety  of  techniques 
used  in  the  field  of  image  reconstruction,  or  tomography.  It  is  used  extensively  in  the 
scientific  and  medical  fields,  allowing  for  an  internal  profile  of  an  object  through 
the  use  of  external  viewing  probes.  The  development  of  the  theoretical  basis  for 
tomographic  methods  is  credited  to  Alan  Cormack  in  the  mid  1960’s  [4,  5].  His 
research,  and  its  subsequent  application  by  Godfrey  Hounsfield  resulted  in  the  creation 
of  computer-assisted  tomography,  or  the  CAT  scan,  and  a  Nobel  Prize  in  Physiology 
and  Medicine  in  1979.  While  the  CAT  scan  method  involves  radiating  an  object 
with  an  external  X-ray  source  and  combining  2-D  profiles  into  a  single  3-D  image, 
the  tomographic  reconstruction  of  plasma  is  slightly  different.  Since  the  plasma  is  a 
radiating  body,  the  measuring  probe  can  be  a  passive  collector  instead  of  an  irradiating 
source.  The  possible  lack  of  symmetry  also  requires  the  use  of  multiple  detectors  in 
various  viewing  axes  to  deconvolve  the  plasma  profile.  Using  multiple  photodiodes 
in  an  array  introduces  the  possibility  of  generating  a  profile  of  the  visible  emissivity 
along  the  viewing  chords.  The  ZaP  flow  Z-pinch  makes  use  of  two  photodiode  arrays 
to  measure  chord-integrated  emissivity,  as  well  as  tomographically  create  one  and  two 
dimensional  radial  profiles  of  the  plasma  interior. 

3.1  Plasma  Bremsstrahlung  and  Line  Radiation 

In  order  to  gain  an  accurate  understanding  of  the  plasma  from  tomographic  data, 
it  is  important  to  determine  exactly  what  sort  of  light  the  photodiodes  detect,  their 
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relationship  to  density  and  energy,  and  how  they  affect  the  resulting  current.  Energy 
radiates  from  the  Z-pinch  in  many  different  forms,  but  for  this  experiment,  the  two  of 
greatest  interest  are  bremsstrahlung  and  line  radiation.  Bremsstrahlung  or  braking 
radiation  results  from  the  interaction  of  ions  and  electrons  undergoing  Coulomb  colli¬ 
sions.  As  fast  moving  electrons  approach  the  relatively  stationary  ions,  the  long-range 
Coulomb  field  deflects  them  into  a  hyperbolic  orbit.  This  deflection  is  determined  from 
the  point  of  closest  approach  in  the  absence  of  any  electromagnetic  forces,  called  the 
impact  parameter  (6),  and  the  deflection  angle  (9),  governed  by  the  following  equation 

[7] 


2 (b/b0 


sin  9  —  — zA_Z_El_  (3.1) 

1  +  ( b/bo ) 

where  bo  is  the  relative  impact  parameter  when  deflection  is  90°.  As  the  electron 
orbits  the  ion  during  Coulomb  collision,  it  accelerates,  resulting  in  the  production  of 
electromagnetic  radiation  called  bremsstrahlung.  The  power  radiated,  W,  from  the 
interaction  between  acceleration  and  the  Coulomb  force  is 


W  = 


2^eV6 


(3,2) 


3  (47re0)3  m2c3r4 

In  this  equation,  Ze g-  is  the  effective  ion  charge,  e  the  electron  charge  in  Coulombs, 
eo  the  permittivity  of  free  space,  m  the  ion  mass,  c  the  speed  of  light,  and  r  the 
instantaneous  radius  of  the  electron  from  the  ion.  Integrating  this  value  over  the 
interaction  time  produces  the  following  amount  of  emitted  radiation  per  collision. 


fbrad 


^eV6 


(3.3) 


3  (47reo)3  m2c3vb3 

where  v  is  the  electron  velocity  relative  to  the  ion  (normally  taken  as  the  electron 
thermal  speed).  The  summation  of  this  equation  for  all  Coulomb  collisions  represents 
the  total  amount  of  energy  emitted  from  the  plasma  from  bremsstrahlung.  Manipu¬ 
lating  this  equation  yields  the  power  radiated  per  unit  volume  of  the  plasma,  given 
by 


r\ 


br 


2n  2nineZ^se6Te/2 
31/2  (47re0)3  m3/2c3h 


(3.4) 


17 


where  Te  is  the  electron  temperature,  h  is  Planck’s  constant,  and  nl  and  ne  are  the  ion 
and  electron  number  densities  respectively.  Therefore,  bremsstrahlung  from  a  plasma 
increases  linearly  with  temperature  and  density  due  to  the  higher  rate  in  Coulomb 
collisions  per  unit  volume. 

Line  radiation  on  the  other  hand,  arises  from  the  transitions  of  bound  electrons 
between  energy  levels.  The  power  produced  per  unit  volume  relates  to  the  combi¬ 
nation  of  electron  density  and  high-Z  ion  density,  as  well  as  electron  temperature 
[7].  The  distinct  energy  levels  of  different  ions  result  in  line  radiation  occurring  at 
specific  wavelengths.  This  allows  for  the  identification  of  various  impurities  inside 
the  vacuum  chamber  of  the  experiment  that  superimpose  over  the  bremsstrahlung 
radiation.  Figure  3.1  is  the  radiation  measured  on  the  ZaP  experiment  for  a  5  kV  ca¬ 
pacitor  bank  charge.  The  greatest  line  radiation  contributions  come  from  ions  such  as 
N-III,  O-II,  O-III,  Ha  (H-I),  B-III,  and  C-III.  With  the  exception  of  Ha  at  656.3  nm, 
most  other  significant  line  radiation  occurs  between  300-500  nm.  Placing  high-pass  or 
band-pass  Liters  in  front  of  the  collimator  openings  of  the  photodiode  array  therefore 
significantly  reduces  the  viewed  line  radiation. 

3.2  Centroid  Calculation 

Given  a  linear  array  of  photodiodes,  it  is  very  simple  to  calculate  the  center  of  ernis- 
sivity  of  the  viewed  radiation  from  the  chord-integrated  data.  The  following  equations 
govern  the  location  of  the  centroid  in  two  dimensions. 

f  x  dArea  f  y  dArea 

x  Area  ’  y  Area 

Since  data  from  the  photodiodes  is  taken  at  discrete  points  along  the  array  axis,  it 
is  easiest  to  determine  the  overall  centroid  by  summing  the  centroid  components  of 
many  smaller  subsections  using  matrix  arithmetic.  In  this  case,  the  centroid  is  the 
summation  of  individual  chord  spacing  A ay  times  the  chord  position  ay  and  signal  Ei 


(3.5) 
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ZaP  Plasma  Emission  Spectra 


Figure  3.1:  ZaP  line  radiation.  B-III  (206.6  nm),  C-III  (229.7  nm),  and  N-III  (409.7 
and  413.4  nm)  represent  the  greatest  contributors  of  line  radiation  in  the  ZaP  plasma. 


divided  by  the  respective  area  under  the  signal  curve. 

c  _  Yh=i  A XjXjEj 

Ylf=i a xiEi 


(3.6) 


The  MATLAB  code  centroid. m  used  for  centroid  calculation  is  included  in  Ap¬ 
pendix  A. 4.  These  centroid  data  are  important  since  they  provide  a  comparison 
between  plasma  radiation  and  current.  Comparing  these  data  to  the  appropriate 
x-  or  ?/axis  component  of  the  center  of  current  data  from  the  magnetic  probes  also 
provides  valuable  insight  into  the  characteristics  of  the  ZaP  Z-pinch. 
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3.3  Centroid  Error  Analysis 

While  the  calculation  of  the  centroid  is  somewhat  straight-forward,  there  are  many 
uncertainties  inherent  to  experimental  analysis.  Standard  deviation  measurements 
indicate  the  degree  of  accuracy  of  the  centroid  calculation.  Standard  deviation  (a) 
is  a  measurement  of  dispersion  of  the  centroid  due  to  error  and  is  governed  by  the 
following  equation  [2] , 


a  = 


(. EiXiAxi 


1 

2 


(3.7) 


Ei,  A Xi,  and  xt  are  the  measured  signal,  spacing,  and  position  of  each  chord.  6Ei 
and  Sxi  represent  the  uncertainty  of  these  signals  and  locations.  In  this  case,  the  well 
known  chord  locations  and  spacings  limit  the  relative  error  contribution  from  these 
terms  to  under  1%.  The  signal  variation  term  however,  contributes  a  higher  error 
of  ~  5%.  The  error  from  this  term  is  especially  apparent  for  the  flat  profiles  seen 
during  formation  and  post-quiescent  periods  in  the  pinch,  where  variations  between 
low  Ei  values  are  much  more  noticeable.  Measurements  of  the  centroid  standard 
deviation  determined  by  the  MATLAB  code  centroiderror  .m  listed  in  Appendix  A. 5 
are  plotted  as  error  bars  in  Figure  3.2. 

Because  centroid  location  plays  an  important  role  in  future  data  analysis,  it  is  im¬ 
portant  to  remove  data,  points  where  there  are  noticeable  differences  between  emission 
and  the  calculated  centroid.  The  centroid. m algorithm  removes  data  points  from  flat 
profiles  and  where  the  calculated  centroid  is  beyond  the  outer  channels  of  the  array. 
If  the  maximum  signal  value  of  the  array  is  less  than  3  times  the  average  of  the  two 
outside  channels,  the  data  are  thrown  out.  This  condition  comes  from  comparing  the 
movement  of  the  emission  profile  to  the  calculated  centroid,  as  shown  in  Figure  3.3. 
Further  results  for  center  and  emissivity  comparisons  are  listed  in  Section  5.3  and  5.4. 
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Figure  3.2:  Centroid  of  ZaP  emissivity  and  standard  deviation.  The  distances  from 
each  centroid  location  indicated  by  the  error  bars  represents  the  standard  deviation 
given  by  Equation  3.7,  giving  upper  and  lower  bounds  to  the  emissivity  location.  The 
horizontal  dashed  lines  represent  the  center  of  each  channel  along  the  PDA  viewing 
axis.  For  this  experiment,  the  standard  deviation  of  the  emissivity  centroid  was 
~  1  channel  spacing  during  quiescence  (shown  between  39.2  and  75.4  /zs  indicated  by 
the  vertical  blue  dashed  lines).  When  the  emission  profile  flattens  out  later  in  time, 
the  standard  deviation  of  the  centroid  increases  to  ~  2  channel  spacings. 


3.4  Chord-Integrated  Emissivity  vs.  Density 

Chord-integrated  emissivity  from  the  photodiodes  is  also  useful  for  comparison  to  the 
Z-pinch  electron  density.  The  radiation  per  unit  volume  for  each  unit  of  frequency  is 


21 


Shot #50316017 


Time  (ps) 


Figure  3.3:  Truncation  criteria  for  centroid  data  due  to  emission  movement  and  flat 
profiles.  The  solid  line  represents  the  complete  centroid  data  for  this  given  emission. 
When  the  profile  at  a  given  time  step  is  not  sufficiently  peaked  or  is  outside  the 
viewing  area  of  the  array,  the  data  are  removed.  Only  the  remaining  points,  shown 
in  red,  are  used  for  comparison  and  further  data  analysis.  The  emission  profile  shown 
here  is  saturated  at  35%  of  the  maximum  value  for  viewing  purposes. 


where  P  is  power  per  electron  and  /  (v)  the  electron  distribution.  If  /  (v)  is  isotropic, 
emissivity  in  terms  of  electron  and  ion  density  defined  by  the  term  j  (u)  is  [13] 


4vr j  (u)  =  neriiZ^s 


32vr2 


2  m 


1/2 


e-hv,Tg. 


(3.9) 


^47re0 /  3\/3 m2c3  \^T / 

This  equation  represents  the  spectral  power  emitted  over  a  full  spherical  solid  angle 
per  unit  of  frequency  per  unit  volume.  ne  and  n*  represent  number  densities,  Zes 


the  effective  ion  charge,  e  the  electron  charge  in  Coulombs,  6q  the  permittivity  of  free 
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space,  c  the  speed  of  light,  T  the  temperature  at  thermal  equilibrium,  hv  the  quantum 
energy,  and  g  the  Maxwell-averaged  Gaunt  factor.  Taking  the  emissivity  over  areas  of 
low  line  radiation  and  low  photon  energy  reduces  this  equation  to  simply  the  number 
densities  multiplied  by  a  constant  so  that  E  oc  n\  [12].  Therefore,  by  viewing  the 
wavelength  region  of  the  plasma  emission  with  only  bremsstrahlung,  it  is  possible  to 
determine  a  general  density  profile  of  the  pinch.  The  best  spectra  determined  from 
Figure  3.1  for  this  condition  fell  between  510-580  urn. 

Unfortunately,  these  data  are  of  little  use  for  a  conclusive  calibration  of  the  inter¬ 
ferometer  using  the  PDA,  as  they  view  different  volumes  of  the  plasma  and  therefore 
cannot  determine  a  correct  scaling  factor.  However,  by  assuming  axisymmetry  of  the 
plasma,  it  is  possible  to  make  general  comparisons  between  the  two  measurements. 
Results  for  density  calculations  are  listed  in  Section  5.2. 

3.5  Radial  Profile  Deconvolution 

One  of  the  more  exciting  uses  of  chord-integrated  data  is  for  tomographic  deconvolu¬ 
tion  of  linear  chord  data  into  a  local  radial  profile.  By  using  one  PDA  with  a  linear 
array  of  photodiodes,  it  is  possible  to  model  the  plasma  as  a  combination  of  indepen¬ 
dent  and  piecewise  constant  cylindrical  shells.  This  is  accomplished  using  the  well 
known  Abel  inversion  technique  [1] . 

3.6  Abel  Inversion  and  One  Dimensional  Theory 

The  emissivity  viewed  by  each  chord  is  merely  the  integration  of  all  the  radial  profiles 
e  (r)  over  the  viewed  length  x,  or  E  (y)  =  f  e  (r)  dx.  Taking  advantage  of  the  geometry 
of  the  path  lengths  through  the  plasma  allows  this  integral  to  reduce  to  only  chord  y 
and  radial  r  geometry, 


E{y)  =  2 


(3.10) 
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Inverting  this  relationship  to  solve  for  e  (r)  requires  determining  the  change  in  chord- 
integrated  emissivity  with  respect  to  y,  resulting  in  the  following  equation, 


e  (r )  = - 

7T 


E'(v) 


v'v 


2  _  r2 


dy. 


(3.11) 


As  there  is  only  one  viewing  array  along  a  single  axis,  it  is  important  to  note  that  an 
Abel  profile  assumes  that  the  plasma  exhibits  circular  symmetry. 

With  multiple  viewing  chords  and  discrete  data,  it  is  easier  to  solve  the  Abel 
inversion  through  the  use  of  linear  algebra  and  a  numerical  methods  program  such  as 
MATLAB.  The  Abel  inversion  follows  the  general  equation  Le  =  E.  where  L  is  a  path 
length  matrix,  e  the  radial  emissivity,  and  E  the  chord-integrated  emissivity.  This 
method  assumes  that  the  radial  profile  is  piecewise  constant  over  each  shell.  Figure  3.4 
shows  a  two  chord  example  with  Equation  3.12  as  its  matrix  representation. 


(3.12) 


In  MATLAB ,  the  inversion  takes  advantage  of  the  built-in  least  squares  solver  to 
increase  computational  speed.  While  the  path  length  matrix  in  the  one-dimensional 
Abel  case  is  square,  and  therefore  exactly  determined,  a  least  squares  data-fitting 
routine  proved  advantageous  in  the  overdetermined  two-dimensional  deconvolution 
discussed  later  in  this  chapter.  In  the  square  matrix  case  however,  this  solver  computes 
the  solution  using  Gaussian  elimination.  A  listing  of  the  MATLAB  routine  abell.m 
is  included  in  Appendix  A. 6. 


3.7  Simulated  Data  and  Error  Analysis  for  1-D  Radial  Profiles 

Discrepancies  are  expected  between  any  numerical  solution  and  the  exact  analytic  so¬ 
lution.  In  an  attempt  to  reach  a  convergent  solution  to  a  linear  problem,  assumptions 
and  approximations  are  always  made  [20].  Geometric  considerations  for  the  photodi¬ 
ode  array,  such  as  chord  spacing,  also  play  a  factor  in  the  accuracy  of  the  results.  In 
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y 


Figure  3.4:  Abel  inversion  path  length  geometry 


order  to  get  a  better  understanding  of  these  influences  on  the  Abel  algorithm,  it  is 
important  to  check  the  deconvolved  data  against  a  known  profile. 

The  MATLAB  Abel  inversion  algorithm  used  to  deconvolve  the  PDA  chord- 
integrated  data  takes  advantage  of  the  centroid  finder  used  in  Section  3.2  to  divide 
the  array  profile  into  left  and  right  sections.  Chord-integrated  data  are  determined 
from  known  radial  profiles,  and  a  path  length  matrix  is  created  from  the  known  chord 
spacing.  These  data  are  then  inverted  and  solved  using  Gaussian  elimination.  Peaked, 
hollow,  centered,  and  offset  profiles  were  all  considered  in  this  analysis.  Table  3.7  is  a 
list  of  several  different  synthetic  profiles  and  their  governing  equations.  These  profiles 
correspond  to  the  inversions  shown  in  Figure  3.5. 

Increasing  the  number  of  chords  viewing  the  plasma  decreases  the  variation  be- 
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(g) 


Radius  (cm) 


Radius  (cm) 


Figure  3.5:  Synthetic  Abel  profiles  for  16  and  32  chords.  These  correspond  to  the 
equations  in  Table  3.7.  The  black  line  represents  the  synthetic  profile,  the  blue  line 
the  chord-integrated  profile,  and  the  red  line  the  returned  Abel  profile.  The  points 
on  these  lines  represent  the  chord  locations. 
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Table  3.1:  Synthetic  profiles  and  governing  equations  for  16  and  32  chord  examples 
used  in  Figure  3.5.  Chords  on  each  end  of  the  profile  are  set  to  zero  in  order  to  force 
the  inverted  profile  to  zero  at  the  edges. 


Profile 

Equation 

Flat  profile 

(Fig.  3.7(a)  and  3.5(b)) 

e  (r)  =  1 

r  <  1 

Linear  profile 

(Fig.  3.5(c)  and  3.5(d)) 

e  (r)  =  1 

—  r 

r  <  1 

Parabolic  profile 

(Fig.  3.5(e)  and  3.5(f)) 

e  (■ r )  =  1 

—  r2 

r  <  1 

y,3 

for  0  <  r  <  0.5 

Hollow  profile 

(Fig.  3.5(g)  and  3.5(h)) 

e  (r)  =  < 

(1  -  r)3 

V 

for  0.5  <  r  <  1.0 

tween  the  initial  and  deconvolved  profiles.  Reducing  the  separation  between  the 
centroid  and  chords  immediately  to  its  right  and  left  also  improved  the  continuity  of 
the  numerical  routine  between  the  two  sides. 


3.8  Centering  Errors  of  the  Abel  Inversion 

While  the  Abel  profiles  shown  in  Figure  3.5  are  all  properly  centered  on  the  given 
symmetric  data,  this  luxury  is  unfortunately  not  always  possible  with  a  moving  Z- 
pinch  plasma.  Although  the  inversion  still  assumes  radial  symmetry,  any  discrepancies 
between  the  center  of  the  inversion  shells  and  the  center  of  the  plasma  profile  can  have 
significant  errors.  For  the  abell.m  inversion  program,  the  inversion  center  is  set  to 
the  centroid  location  determined  by  the  program  centroid. m.  Figure  3.6  shows  the 
discontinuity  at  the  center  of  the  Abel  profile  caused  by  offsetting  the  shell  locations 
by  1  chord  spacing.  The  profile  is  given  by  the  equation  e  (r)  =  1  —  r2,  also  shown 
correctly  centered  in  Figure  3.5(f). 
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Figure  3.6:  A  32  chord,  offset  synthetic  profile  given  by  the  equation  e  (r)  =  1  —  r2. 
The  black  line  is  the  synthetic  radial  profile  centered  on  the  black  dashed  line.  The 
blue  line  is  the  chord-integrated  data  taken  from  the  radial  profile.  The  red  line 
indicates  the  radial  Abel  inversion  profile  with  the  shell  centers  set  at  0  cm  on  the 
x-axis. 


3.9  Two  Dimensional  Tomographic  Deconvolution 

While  the  Abel  inversion  technique  is  useful  for  resolving  profiles  with  radial  sym¬ 
metry,  it  cannot  accurately  model  those  with  azimuthal  variations.  Inverting  a  non- 
symmetric  profile  can  lead  to  discontinuities  between  the  left  and  right  sides.  Rotating 
the  viewing  array  into  different  axes  also  reveals  changing  radial  profiles.  Inverting 
chord-integrated  data  at  cross-sections  where  the  profile  is  not  symmetric  also  leads 
to  discontinuities  between  the  left  and  right  side  sections,  as  shown  in  Figure  3.7. 

Simultaneously  viewing  the  plasma  in  different  orientations  using  multiple  arrays 
allow  for  azimuthal  variations  to  be  considered  in  the  deconvolution  profile.  This 
allows  for  better  resolution  of  the  emissivity  of  the  plasma  relating  to  the  mode 
instabilities  listed  in  Section  2.2.  Fitting  the  chord-integrated  data  of  two  photodiode 
arrays  to  superimposed  Fourier  component  profiles  and  solving  a  similar  numerical 
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(a)  Azimuthally  varying  hollow  profile 


(b)  Y-axis  cross-section 


(c)  X-axis  cross-section 


Figure  3.7:  (a)  Abel  inverted  cross-sections  of  an  azimuthally  varying  profile,  (b)  is 
the  inversion  taken  along  the  x-axis,  and  (c)  along  the  y- axis,  where  the  black  line 
is  the  synthetic  profile,  the  blue  line  the  chord- integrated  data,  and  the  red  line  the 
inverted  profile.  The  dots  represent  the  discrete  data  points.  Note  the  discontinuity  in 
(c)  resulting  from  the  Abel  algorithm  attempting  to  follow  the  non-symmetric  chordal 
data  across  the  left  and  right  sections.  The  hollow  profile  in  this  example  follows  the 

,  .  ,  .  .  .  Jr3  for  0  <  r  <  0.5 

equation  e  (r)  =  en  (r)  +  0.3en  in  cosd,  where  en  (r)  =  < 

V;  J  1  (1  —  r)3  for  0.5  <  r  <  1.0 
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linear  algebra  routine  as  the  one-dimensional  case  allows  for  the  generation  of  a  shell 
profile  of  the  Z-pinch  varying  both  radially  and  azimuthally. 

3.10  Fourier  Component  Theory 

Solving  a  two  dimensional  plasma  profile  involves  using  a  Fourier  series  expansion 
method  developed  by  Cormack  [4,  5].  A  radially  and  azimuthally  varying  plasma 
profile  is  represented  in  polar  coordinates  by  g  (r,  9).  For  a  chord  viewing  a  circular 
plasma  through  a  length  L ,  the  line  integral  of  the  profile  is  simply 

f  ip  A)  =  [  9  (r ,  9)  ds.  (3.13) 

J  L(p,<j>) 

p  represents  the  perpendicular  distance  from  the  geometric  origin  of  the  profile  to  the 
line  L ,  and  0  is  the  angle  from  the  positive  x-axis  to  the  line  measuring  p.  Because 
of  the  viewing  geometry  for  the  ZaP  experiment  and  the  restricted  use  of  only  two 
PDAs,  the  value  of  p  is  equal  to  each  array’s  respective  chord  location  along  the 
viewing  axis.  0  is  also  limited  to  a  value  of  0°  or  90°,  again  depending  on  the  viewing 
axis.  In  addition,  if  g  is  assumed  to  be  piecewise  constant,  this  expands  in  Fourier 
series  to 

l 

9M)=5>(r”>)e""  (3.14) 

n= 0 

where  gn  is  the  piecewise  constant  profile  for  shell  m.  Without  e  ,  this  is  simply 
the  Abel  profile  described  in  Section  3.6.  Expanding  the  complex  em 6  into  sine  and 
cosine  components  defined  by  ein6  =  cos  n9  +  %  sin  n9  results  in 

g  (r,  9)  =  g0  (rm)  +  gi  (rm)  [cos6»  +  z  sin  6>]  .  (3.15) 

This  equation  represents  a  non-symmetric  radial  shell  profile  by  combining  the  de¬ 
convolved  Abel  profile  with  one  multiplied  by  cosine  and  sine  functions. 

As  with  the  Abel  inversion  algorithm,  the  two-dimensional  profile  lends  itself  to 
matrix  computation  and  linear  algebra.  Rewriting  equation  3.15  so  that  each  shell 
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profile  e  (r)  =  eo  +  e\  cos  9  +  e 2  sin0,  the  relationship  between  chord-integrated  data 
and  radial  data  takes  the  form  of  E  =  Me,  where  E  are  chord-integrated  values, 
M  are  the  combined  path  length  and  sine  and  cosine  components,  and  e  are  the  eo, 
ei,  and  e2  values  for  each  shell.  A  criterion  that  at  least  three  chords  pass  through 
the  outer  shell  is  used  to  determine  the  maximum  number  of  shells  deconvolved  for 
a  given  profile.  Smaller  shells  are  added  inside  of  this  one  at  a  distance  equal  to  the 
average  chord  spacing.  This  ensures  that  the  number  of  chords  provided  by  the  PDAs 
will  always  be  greater  than  or  equal  to  the  unknown  eo,  ei,  and  e2  profiles  for  all  the 
shells,  or  rather  that  it  is  either  exactly  or  overdetermined.  An  elementary  one  shell 
example  is  given  in  Figure  3.8,  along  with  the  corresponding  matrix  represenation  in 
Equation  3.16,  where  9  is  the  angle  from  the  positive  rr-axis  to  the  midpoint  of  each 
pathlength  element. 


E11 

Ei  2 

E21 

E22 

2 -bn 
2Ti2 
2L2i 
2-^22 


L\\  cos  0n  +  L\\  cos  (— 0n) 

L\2  COS  012  +  L12  COS  (-012) 
L2 1  COS  (-022)  +  L2\  COS  (  021 ) 
L22  cos  022  +  L22  COS  021 


Ln  sin  0n  +  Ln  sin  (-0n) 
L12  sin  0i2  +  L12  sin  (-0i2) 
L2i  sin  (-022)  +  L21  sin  (-02i) 
L22  sin  022  +  L22  sin  02i 


(3.16) 


Unlike  the  1-D  Abel  case,  this  method  solves  for  all  the  chords  from  both  the  top 
and  side  viewing  arrays  in  one  single  path-length  matrix  instead  of  dividing  them  by 
their  respective  centroid  locations.  Because  the  path-length  matrix  M  can  be  over¬ 
determined  and  therefore  rectangular,  the  solving  algorithm  requires  a  least  squares 
fit  routine  to  calculate  the  best  inverted  profile.  The  MATLAB  routines  TDerror.m 
and  bestf  it  .m  make  use  of  QR  decomposition  with  pivoting  to  converge  to  the  best 
solution.  9  values  in  this  matrix  represent  the  angle  from  the  positive  x-axis  of  the  first 
PDA  to  the  midpoint  of  the  respective  chord  through  an  individual  shell.  This  small 
angle  approximation  becomes  less  pronounced  with  an  increasing  number  of  chords.  A 
listing  of  these  two  MATLAB  routines  and  the  overall  Fourier  deconvolution  algorithm 
FULL2d.m  are  given  in  Appendix  B. 
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Figure  3.8:  2-D  Fourier  deconvolution  path  length  and  angle  geometry.  Line  segment 
midpoints  are  used  to  determine  6  angles  for  cosine  and  sine  components. 


3.11  Simulated  Data  and  Error  Analysis  for  2-D  Radial  Profiles 

Comparisons  of  the  deconvolved  data  to  synthetic  radial  profiles  are  just  as  important 
in  understanding  the  limitations  of  the  2-D  Fourier  deconvolution  method  as  they  were 
for  the  1-D  Abel  inversion  algorithm.  Validating  the  code  involved  using  data  from 
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azimuthally  varying  profiles,  both  centered  and  offset  with  respect  to  the  photodiode 
arrays.  Results  exposed  several  different  factors  that  play  an  important  role  in  the 
amount  of  error  generated  in  the  2-D  numerical  method. 

In  order  to  minimize  the  approximations  and  converge  to  a  solution,  the  least 
squares  fit  solver  attempts  to  find  the  minimum  norm,  || a;  ||  2,  with  no  rank  deficiency. 
A  rank  deficient  matrix  indicates  that  several  channels  are  not  linearly  independent, 
meaning  that  there  are  not  enough  data  points  to  resolve  unknowns  in  some  of  the 
shell  profile.  This  occurs  primarily  from  not  having  enough  independent  chords  for 
either  the  outer  or  inner  shell  profiles.  Minimizing  these  conditions  involved  setting 
the  maximum  shell  radius  to  the  criteria  set  in  Section  3.10  and  by  guaranteeing  a 
minimum  of  four  chords  traverse  the  inner  shell.  Although  only  three  chords  are 
required  to  solve  for  the  three  unknown  eo,i,2  values,  when  a  chord  passes  very  near 
the  origin  of  the  inner  shell,  the  larger  jump  in  path-length  sizes  and  the  proximity  of 
their  9  angles  to  0  or  90°  can  skew  the  data.  Increasing  the  inner  shell  radius  to  allow 
an  extra  chord  to  pass  through  reduces  the  chances  that  this  condition  will  exist. 

Correct  geometric  centering  of  the  inverted  profile  on  the  plasma  center  is  also 
a  serious  concern  in  error  analysis.  Since  path-lengths  passing  through  the  origin 
of  the  radial  model  can  distort  the  profile  of  the  inner  shell,  the  algorithm  forces  a 
minimum  perpendicular  separation  between  any  given  chord  and  the  chosen  center. 
Unfortunately,  since  the  entire  shell  model  geometry  revolves  around  this  point,  any 
offset  in  the  routine  center  can  also  have  disastrous  consequences  and  result  in  a 
deconvolved  profile  very  different  from  the  synthetic  contour.  A  minimum  error  finder 
determines  the  best  profile  center  by  comparing  the  given  and  the  calculated  chord- 
integrated  data  and  hireling  the  minimum  y2  location,  y2  is  a  measurement  of  the 
goodness  of  fit  of  least  squares  data,  and  is  governed  by  the  following  equation  [2] , 

X2  =  XI  (“2  [Vi  ~  V  Ou)]2)  •  (3.17) 

The  standard  deviation  (a)  squared  term  is  the  variance,  or  expected  error  in  the 
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Table  3.2:  Synthetic  profiles  for  2-D  Fourier  analysis  used  in  Figures  3.9(f) — 3.11(f). 
An  azimuthal  variation  is  added  to  e  (r)  by  sine  and  cosine  components  30%  the  size 
of  the  axisymmetric  data. 


Profile 

e  (r) 

=  e  +  0.3e  cos  9  +  0.3e  sin  9 

Parabolic  (Fig.  3.9(a),  3.10(a),  3.11(a)) 

Hollow  (Fig.  3.9(d),  3.10(d),  3.11(d)) 

e  =  1 

e  =  < 

—  r2  r  <  1 

/ 

r3  for  0  <  r  <  0.5 

(1  -  r)3  for  0.5  <  r  <  1.0 

signal  of  each  channel.  yl  represents  the  measured  signal  and  y  (ay)  is  the  returned 
calculated  signal.  When  the  fitting  routine  is  a  good  approximation  of  the  given  syn¬ 
thetic  data,  the  value  of  y2  is  on  the  order  of  1.  Contour  plots  of  several  azimuthally 
varying  profiles  and  their  respective  y2  data  are  shown  in  Figure  3.9.  Table  3.11  is 
a  list  of  the  equations  used  to  generate  the  synthetic  profiles.  These  data  show  that 
the  width  of  the  contour  of  y2  =  1  is  always  much  less  than  the  spacing  between 
the  viewing  chords,  meaning  that  the  final  position  of  the  profile  center  chosen  by 
the  algorithm  is  very  well  known  compared  to  the  chord  spacings.  The  MATLAB 
algorithm  FULL2d.m  and  subroutines  TDerror  .m  and  bestf  it  .m  calculate  the  decon¬ 
volved  2-D  profile  at  the  minimum  y2  position.  Keeping  track  of  location  over  time 
also  allows  the  algorithm  to  track  the  plasma  profile  movement  off  the  2-axis,  as  shown 
in  Figure  3.10. 

The  TDerror  .m  and  bestf  it .  m  routines  also  compensate  for  cases  where  the  large 
variation  in  sine  and  cosine  contributions  lead  to  incorrect  profiles.  This  involves 
constraining  the  value  of  the  axisymmetric  eo  term  by  the  e.\  and  e 2  data  according 
to  the  following  equation, 

|eo|  >  \Je\  +  e22.  (3.18) 

Without  this  constraint,  the  sine  and  cosine  components  of  the  deconvolution  can 
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cause  the  routine  to  return  a  negative  emissivity  profile.  These  eo,  ei,  and  e2  values 
are  used  to  determine  the  calculated  e  (r)  values.  Multiplying  e  (r)  and  the  path- 
length  matrix  M  produces  the  calculated  signals  y  (xi)  used  in  y 2  minimization. 

In  addition,  random  noise  is  added  to  the  algorithm  to  test  the  effects  of  digitizer 
noise  in  profile  reconstruction.  This  is  accomplished  by  adding  random  noise  with  a 
known  magnitude  to  the  chord-integrated  data.  The  at  term  in  Equation  3.17  is  set 
to  this  magnitude  in  order  to  match  the  known  error.  During  experimentation,  a,  is 
determined  by  calculating  the  bit  noise  of  each  digitizer  channel.  Figure  3.11  shows 
the  effects  of  adding  a  5%  error  to  the  input  signal. 
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Figure  3.10:  2-D  synthetic  and  deconvolved  profiles  with  x  and  y  =  -0.5  cm  offset.  The  algorithm  for  the  circular 
profile  set  the  reconstruction  center  at  x  and  y  =  -0.43  cm.  The  deconvolved  hollow  profile  is  centered  at  x  and  y  = 
-0.53  cm.  Both  profiles  are  still  within  a  single  chord  spacing  of  the  synthetic  data. 
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Figure  3.11:  2-D  synthetic  and  deconvolved  profiles  with  added  noise  equivalent  to  5%  of  the  maximum  signal.  The 
noise  has  an  average  of  zero  and  a  standard  deviation  of  1.  Again,  the  goodness  of  fit  determined  by  y2  =  1  varies  by 
less  than  the  chord  spacing. 
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Chapter  4 

PHOTODIODE  ARRAY  DESIGN  AND  CALIBRATION 

Using  a  photodiode  to  detect  plasma  emissivity  is  a  straightforward  process  that 
simply  involves  pointing  it  at  the  light  source  and  measuring  the  resulting  current. 
Unfortunately,  taking  useful  tomographic  data  requires  a  great  deal  more  design  con¬ 
sideration  and  calibration.  An  array  of  photodiodes  collects  light  from  different  re¬ 
gions  of  the  plasma.  Collimation  of  these  photodiodes  reduces  the  visible  solid  angle 
and  provides  meaningful  data  based  upon  the  array  configuration.  This  requires  a 
well  known  array  geometry  and  defined  viewing  angles  that  remove  reflected  light. 
The  vacuum  tank  window  design  and  inner  electrode  openings  further  restrict  the 
area  of  the  measured  plasma.  Components  of  the  overall  photodiode  array  used  in 
experimentation  are  given  in  Figure  4.1.  Knowledge  of  the  spectral  response  of  the 
photodiode  chip,  individual  signal  amplification  and  attenuation,  frequency  response 
and  noise  filtering  are  all  required  to  accurately  interpret  the  resulting  data. 

4.1  Photodiode  Chip  Design,  Spectral  Response,  and  Biasing  Board 

As  previously  mentioned  in  Section  2.1,  an  aligned  array  of  photodiodes  allows  for 
the  collection  and  generation  of  an  emissivity  profile  of  the  Z-pinch  plasma.  The 
ZaP  experiment  makes  use  of  four,  16-element  silicon  photodiode  arrays  made  by  the 
Hamamatsu  Corporation.  The  S5668-02  series  proved  ideal  in  providing  a  quick  signal 
rise  time  and  linear  array  configuration  for  tomographic  analysis  over  the  visible  light 
spectrum.  Each  PDA  uses  two  aligned  chips  to  create  a  32  channel  linear  array.  Since 
the  photodiodes  all  have  a  common  cathode,  they  produce  a  small  negative  closed  loop 
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TANK  CONNECTION  —  gg5>4.62  in  diameter 
FLANGE  o.20  in  thickness 

0.332  in  bolt  holes 
2.50  x  0.75  in  opening 


Figure  4.1:  Components  of  the  photodiode  array 
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current  when  exposed  to  a  light  source.  Each  active  element  has  a  1.175  mm  width 
and  2  mm  height  along  the  viewing  axis.  The  mounting  PC  board  is  25  mm  (1  in) 
in  length,  allowing  for  easy  mounting  and  alignment.  A  schematic  of  the  photodiode 
chip  geometry  is  shown  in  Figure  4.2  [10].  Table  4.1  lists  typical  electrical  and  optical 
characteristics  of  the  photodiode  chip. 


NOTE  DIMS  IN  Ml 

Cowion  Cathode  Pin 


1,0 


■FeNiCo  PIN 


•  Glass  Epoxy  Board 
G-10  < Black) 

■Epoxy  Resin 


3,5 


Figure  4.2:  Hamamatsu  S5668-02  silicon  photodiode  chip.  Dimensions  are  in  milli¬ 
meters. 


As  shown  in  Table  4.1,  the  S5668-02  series  photodiodes  have  wavelength  detection 
limits  between  320-1100  nm,  with  a  peak  sensitivity  of  960  nm.  The  photocurrent 
produced  from  incident  light  between  these  boundaries  varies  according  to  wavelength 
along  a  spectral  response  curve,  shown  in  Figure  4.3  [10].  The  higher  order  lines  seen 
in  the  spectra  plots  represent  the  light  scattered  in  frequency  space  by  the  diffraction 
grating  of  the  ICCD  spectrometer.  They  do  not  contribute  to  the  actual  amount 
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Tabic  4.1:  Important  electrical  and  optical  characteristics  for  the  Hamamatsu  S5668- 
02  silicon  photodiode  chip. 


Parameter 

Symbol 

Unit 

Spectral  Response  Range 

A 

320-1100 

nm 

Peak  Sensitivity  Wavelength 

A  p 

960 

nm 

Photo  Sensitivity 

S 

See  Figure  4.3(a)  for  curve 

A/W 

Rise  Time 

tr 

0.1 

ps 

of  line  radiation  in  the  plasma  and  can  therefore  be  ignored.  As  previously  seen 
in  Figure  3.1,  B-III,  C-III,  and  N-III  ions  produce  a  great  deal  of  line  radiation. 
Fitting  the  response  curve  to  the  plasma  radiation  plot  reveals  that  B-III  and  C- 
III  are  completely  filtered  out  by  the  photodiode  chips.  N-III  remains  a  significant 
contributor  to  cmissivity  for  these  pulses,  however  higher  wavelength  line  radiation 
now  make  up  a  larger  proportion  of  the  light  detected. 

Since  plasma  instabilities  in  the  ZaP  sheared  flow  Z-pinch  occur  on  the  order  of 
microseconds  or  faster,  it  is  imperative  to  increase  the  response  time  of  the  photodiode 
circuit  to  this  same  order  of  magnitude.  A  biasing  circuit  running  off  of  a  9  V  battery 
accomplishes  this  by  introducing  a  reverse  bias  current  on  the  photodiodes.  This 
current  replenishes  the  electron  hole  pairs  in  the  silicon  crystal  depletion  layer  between 
the  p  and  n  sections  of  the  photodiode,  thereby  reducing  the  signal  current  delay.  This 
board  mounts  directly  onto  the  back  of  the  photodiode  chips,  as  shown  in  Figure  4.4. 


4.2  Collimator  Design 

Alignment  and  collimation  of  the  photodiode  chips  is  accomplished  using  a  DuPont 
plastic  polymer  called  Delrin.  Delrin  proved  advantageous  because  of  its  durability, 
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ZaP  Plasma  Emission  Spectra 


Wavelength  X  (nm) 

(a) 


ZaP  Plasma  Emission  Spectra 
as  seen  by  the  Photodiodes 


Figure  4.3:  Effects  of  photodiode  sensitivity  on  viewed  emission  spectra,  (a)  is  the  full 
emission  spectra  with  the  the  photodiode  crystal  photo  sensitivity  curve  overplotted, 
(b)  is  the  ZaP  spectra  as  seen  by  the  photodiodes  due  to  the  response  curve.  N-III 
remains  the  top  supplier  of  line  radiation.  Higher  order  lines  seen  in  the  spectrometry 
measurement  may  be  ignored. 
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Figure  4.4:  PDA  biasing  board  mounted  on  the  photodiode  chips.  Each  set  of  twisted 
pair  wire  connects  to  one  of  the  thirty-two  channels  from  the  photodiodes. 


ease  to  machine,  smooth  surfaces,  and  its  availability  in  a  black  opaque  color  that 
would  minimize  reflected  light.  Four  pieces  of  Delrin  plastic  were  used  in  the  con¬ 
struction  of  each  of  the  two  collimators,  and  most  of  the  geometric  restrictions  were 
a  direct  result  of  limits  from  the  window  mounting  on  the  ZaP  vacuum  tank  and  2  in 
diameter  outer  electrode  view  ports.  A  previous  sixteen  element  photodiode  array 
device  used  on  the  ZaP  experiment  offered  a  template  for  the  base  design  of  the  cur¬ 
rent  arrays.  An  exploded  view  of  the  current  collimator  design  is  given  in  Figure  4.5 
and  a  color  enhanced  schematic  of  the  combined  assembly  is  shown  in  Figure  4.6. 

The  photodiode  chips  are  housed  in  a  recess  in  a  3.6  x  1.75  x  0.25  in  Delrin  cradle 
piece.  The  sides  of  the  Delrin  opening  provide  horizontal  alignment  of  the  chips 
above  thirty-two  0.032  in  diameter  holes.  Each  hole  is  centered  0.062  in  apart  along 
the  horizontal  axis  with  tolerance  levels  minimized  to  approximately  0.001  in  through 
the  use  of  a  digital  display  endmill.  The  chips  are  secured  by  backing  socket  pins  and 
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8-32  Collimator 
screw  hole 


RUBBER  SPACER 


8-32  bolthole 


screw  hole 


SMOOTH  SIDE  PIECE 


TOP  PLATE 


8-32  Biasing  Board 
-  screw  hole 


Alignment  Pin 

0.187  In) 


SLOT  PIECE 


Collimator  Slots 
(width  0.040  In) 


Aiignme 
holes  tala.  0.187  In) 

8-32  Collimator 
screw  hole 


Collimator  Holes 

(dla.  0.032  In) 


Figure  4.5:  Unassembled  photodiode  array  collimator 
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8-32  Collimator 
connection  screw 


8-32  L-flange  8-32  Collimator 
screw  hole  connection  screw 


Figure  4.6:  Assembled  and  aligned  photodiode  array  collimator 


3.28  in 
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the  biasing  board,  which  attaches  to  the  Delrin  cradle  using  four  8-32  screws.  The 
cradle  assembly  sits  atop  two  interconnecting  side  pieces,  a  large  smooth  side  piece 
and  the  slot  piece.  The  slot  piece  has  thirty-two  0.040  in  wide  slots  with  0.022  in 
ridges  and  is  3.6  x  0.91  x  2.78  in  in  size.  Alignment  between  the  cradle  holes  and  the 
slot  piece  occurs  through  the  use  of  two  0.187  in  alignment  pins,  secured  in  specially 
drilled  holes.  Four  8-32  bolts  connect  the  sides  together  like  puzzle  pieces.  A  small 
0.020  in  thick  opaque  rubber  piece  sits  between  the  smooth  side  and  the  slots.  The 
compression  force  created  by  the  bolts  ensures  that  this  rubber  piece  will  close  any 
gaps  between  neighboring  slots,  reducing  the  chance  that  light  will  “leak”  from  one 
channel  to  another.  The  side  pieces  combine  to  produce  the  same  length  and  width 
profile  as  the  cradle  piece.  A  “top  plate”  with  the  same  dimensions,  alignment  pins, 
and  hole  configuration  as  the  cradle  attaches  to  the  other  end  of  the  side  pieces 
to  create  a  rectangular  box  collimator.  A  smaller  recess  was  cut  into  this  piece  to 
allow  for  easier  drilling  of  the  collimation  holes,  and  to  minimize  difference  in  pitch 
caused  by  the  drill  press.  In  all,  the  four  pieces  combine  to  offer  approximately  2.9  in 
of  collimator  length  from  the  top  plate  opening  to  the  photodiode  crystal.  Overall 
design  tolerance  for  the  assembly  was  approximately  0.003  in. 


Since  the  signal  detected  by  the  PDA  is  directly  proportional  to  the  amount  of 
light  striking  each  photodiode  crystal,  it  is  important  to  know  the  viewing  solid  angle 
of  each  channel.  Although  the  collimator  design  highly  parallelizes  the  light,  by  the 
time  the  light  cones  reach  the  center  of  the  tank  (7  in  from  the  vacuum  chamber 
window),  the  viewing  area  created  by  the  0.032  in  (0.08  cm)  diameter  collimator  hole 
balloons  to  a  diameter  of  0.19  in  (0.49  cm).  As  shown  in  Figure  4.7,  this  creates  a 
total  viewing  area  of  2.11  x  0.19  in  (5.37  x  0.49  cm)  where  adjacent  channels  overlap. 
This  overlap  acts  to  smooth  the  chord-integrated  profile  recorded  by  the  PDA.  Using 
fewer  channels  on  the  PDA  allows  for  the  generation  of  a  more  discrete  profile. 
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Figure  4.7:  Collimator  viewing  areas  of  the  ZaP  plasma.  All  dimensions  are  given  in 
centimeters,  (a)  shows  the  collimator  hole  openings  and  (b-d)  the  channel  viewing 
angles  of  a  1  cm  radius  plasma  7  in  from  the  collimator  circular  connection  flange, 
(b)  shows  8  equally  spaced  channels,  (c)  16  channels,  and  (d)  all  32  channels. 


4.3  Amplifier  Board  Design 

While  the  collimator  excels  at  highly  parallelizing  the  viewed  light,  it  has  the  unfor¬ 
tunate  consequence  of  significantly  reducing  the  photodiode  output  current.  Back¬ 
ground  noise  is  therefore  much  more  significant,  and  special  care  must  be  taken  to 
shield  collection  circuits  from  the  lab  environment  and  filter  out  unwanted  noise.  The 
amplifiers  also  match  the  output  impedance  of  the  photodiodes  to  the  50  fl  input  of 
the  digitizers. 

Each  photodiode  in  the  array  has  an  amplification  circuit  consisting  of  two  op- 
amps;  a  transimpedance  amplifier  and  a  signal  buffer.  Figure  4.8  shows  the  electrical 
schematic  for  a  single  amplifier.  The  AD797  op-amp  provides  signal  amplification 
with  very  little  noise  distortion.  Unfortunately,  it  also  has  a  small  50  mA  peak  output 
current.  A  subsequent  LM7171  buffer  increases  this  output  current  to  a  usable  level 
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Trimmer  Capacitor 


1.5-10p 


Figure  4.8:  Photodiode  amplifier  circuit  diagram 


and  provides  voltage  amplification  of  the  final  signal.  The  47  resistor  following  the 
LM7171  prevents  transmission  reflections  from  the  cables  in  the  lab  and  protects  the 
amplifier  in  the  event  of  an  electrical  short.  The  only  variable  component  of  the  design 
is  a  1.5-10  pF  trimmer  capacitor,  which  allows  for  frequency  filtering  by  correcting 
underdamped  or  overdamped  signal  outputs. 

High  background  noise  levels  in  the  lab  caused  by  the  main  bank  power  supply 
and  plasma  current  interruption  during  post-quiescence  are  of  great  concern  to  the 
photodiode  array  design  as  they  can  interfere  with  the  collection  of  any  meaningful 
data.  The  amplifiers  filter  out  this  noise  by  doubling  as  low  pass  frequency  filters. 
Frequency  response  tests  on  the  amplifier,  as  shown  in  Figure  4.9,  indicate  that  they 
successfully  remove  high  frequency  noise  from  data  with  a  frequency  below  2  MHz 
without  negative  gain.  Because  the  frequencies  of  interest  for  the  ZaP  experiment  are 
well  below  this  2  MHz  corner,  any  remaining  high  frequency  noise  can  be  neglected  by 
smoothing  the  data.  The  amplifiers,  collimator  cradle,  biasing  board,  and  photodiode 


49 


PDA  Frequency  Response 


Figure  4.9:  Photodiode  amplifier  frequency  response  adjusted  to  DC  gain  =  1.  The 
black  points  represent  signal  gain  at  given  input  frequencies.  The  red  line  is  a  rough 
trendline  of  these  points.  These  data  show  that  the  amplifiers  can  filter  out  high 
frequency  noise  without  signal  degradation  below  an  input  frequency  of  2  MHz. 


chips  are  all  shielded  from  the  plasma  magnetic  field  by  a  thin  aluminum  box,  as  shown 
in  Figure  4.10.  The  output  signals  for  the  thirty-two  channels  are  divided  between 
three  DB-25  connectors  on  a  removable  side  panel  of  the  box.  In  order  to  fit  everything 
inside  the  box,  the  amplifier  circuits  are  mounted  onto  PC  boards  and  stacked  using 
8-32  threaded  rod  and  insulating  spacers.  The  sixteen  odd  channels  and  sixteen  even 
channels  are  divided  into  two  groups  of  four  boards,  each  receiving  a  separate  ±12  V 
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DB-25  Connectors 


Amplifier  Power 


Trimmer  Capacitor 
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5  v  -  _ _ —  » 


LM7171  *  •  :W 
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Cbllimator,  Biasing  board,  and  Photodiode  chips 


Figure  4.10:  Stacked  amplifier  boards  inside  the  shielding  aluminum  box 


power  supply  required  by  the  op-amps.  This  division  of  the  input  power  is  important 
because  of  the  high  supply  current  required  by  the  amplifier  boards.  Each  side  of  the 
AD797  and  LM7171  chips  have  peak  input  currents  of  10.5  and  9.5  mA  respectively. 
The  thirty-two  channels  add  up  to  an  overall  peak  supply  current  requirement  of 
1.28  A.  Dividing  the  input  current  into  two  equal  parallel  supplies  reduces  the  risk 
of  overheating  any  of  the  circuit  components.  Thermocouple  tests  conducted  on  the 
PDA’s  during  data  collection  showed  that  element  temperature  levels  remained  well 
below  the  85°  C  (185°  F)  maximum  operating  temperature  using  this  power  supply 
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configuration.  The  external  power  is  also  divided  between  the  DB-25  connectors 
and  twisted  pair  cable  running  out  of  the  laboratory.  A  picture  of  the  overall  PDA 
assembly  mounted  on  the  ZaP  Experiement  is  shown  in  Figure  4.11. 


4.4  Photodiode  Array  Channel  Calibration 

Due  to  variations  in  the  photodiode  chip,  collimator,  biasing  and  amplifier  boards, 
calibration  factors  are  required  to  correlate  the  channels  on  both  arrays  to  one  relative 
standard.  This  is  accomplished  by  relating  the  maximum  signals  detected  by  each 
channel  from  a  constant  light  source.  During  the  testing  and  calibration  phase,  the 
photodiode  array  was  secured  on  a  translational  optics  table  such  that  the  moving  axis 
was  aligned  with  the  photodiode  viewing  axis.  A  Strobotac®  1531-AB  series  strobe 
light  made  by  the  General  Radio  Company  with  a  flash  frequency  of  approximately 
167  Hz  was  placed  behind  a  frosted  glass  window  and  an  opaque  cardboard  covering 
with  a  narrow  1.0  in  slit  window.  The  frosted  glass  and  slit  window  ensured  that  only 
a  small  band  of  uniform  light  was  visible  to  the  PDA.  The  strobe  light  was  placed 
on  the  translational  table  approximately  7  in  from  the  collimator  opening,  roughly 
the  same  distance  it  would  be  from  the  center  of  the  plasma  when  in  use  on  the 
experiment.  Figure  4.12  shows  a  picture  of  the  channel  calibration  setup.  An  external 
±12  V  supply  was  applied  to  power  the  PDA  amplifier  boards.  An  oscilloscope  was 
attached  to  the  DB-25  connectors  to  measure  channel  output  signals  as  the  strobe 
light  translates  along  the  optics  table.  The  peak  to  peak  output  signal  is  recorded  for 
each  channel  when  the  slit  window  was  aligned  directly  over  that  particular  collimator 
hole.  The  movement  of  the  table  also  provides  a  rough  measurement  of  the  distances 
between  the  viewing  angle  centers  of  adjacent  channels.  The  average  distance  of 
0.15-0.16  cm  measured  corresponds  to  the  distances  between  ray  traces  in  Figure  4.7. 
Figure  4.13  shows  an  example  of  the  measured  strobe  light  signal  and  translation 
as  seen  by  the  oscilloscope.  This  testing  also  provided  a  means  to  check  amplifier 
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Figure  4.11:  Photodiode  array  mounted  on  the  ZaP  experiment 
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Amplifier  box 
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Figure  4.12:  Channel  calibration  setup 


(a) 


(b) 


Figure  4.13:  Channel  calibration  measurements  as  seen  by  the  oscilliscope.  Four 
output  channels  are  displayed  in  each  plot.  In  Figure  (a),  the  strobe  light  slit  is 
aligned  over  the  collimator  hole  for  ch.  1.  In  Figure  (b),  the  strobe  light  is  translated 
so  that  the  collimator  hole  for  ch.  4  is  now  aligned  with  the  slit.  The  peak-to-peak 
voltage  measured  for  each  channel  is  used  to  create  calibration  factors. 
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damping  levels  and  circuit  connection  problems  inside  the  amplifier  box.  Over  and 
underdamped  signals  were  adjusted  to  corrected  levels  using  the  trimmer  capacitors 
on  the  amplifier  boards,  and  any  noise  spikes  or  loose  connections  were  investigated 
and  fixed  before  recalibrating  that  channel.  Figure  4.14  illustrates  several  different 
problems  encountered  during  channel  calibration. 


(b) 


(c) 


Figure  4.14:  Amplifier  board  errors.  Figures  (a)  and  (b)  show  a  flatline  and  internal 
noise  signal  respectively.  These  are  indicative  of  a  loose  solder  connection  on  that 
particular  channel.  Figure  (c)  shows  an  underdamped  amplifier  response.  Adjust¬ 
ing  the  trimmer  capacitor  shown  in  Figure  4.10  fixes  this  problem,  however  a  new 
calibration  factor  must  then  be  taken  for  that  channel. 
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The  biasing  and  amplifier  boards  were  also  tested  independent  of  the  collimator 
for  cross  talk  and  frequency  response  characteristics.  With  the  socket  removed  from 
the  photodiode  chips,  a  square  signal  from  a  waveform  generator  was  supplied  to 
the  circuit  as  shown  in  Figure  4.15.  The  frequency  response  of  the  amplifier  circuits 
calculated  from  this  test  is  approximately  8  MHz,  well  above  the  frequency  of  interest 
for  this  experiment.  A  second  oscilloscope  channel  measured  the  effects  of  cross  talk 
between  adjacent  channels  on  the  amplifier  boards.  The  noise  measured  by  the  second 
channel  was  under  1%  of  the  input  square  waveform  and  was  therefore  negligible. 


Figure  4.15:  Amplifier  board  frequency  response  and  channel  cross  talk  measurements. 
Channel  1  measures  a  rise  time  of  127.4  ns  associated  with  the  square  waveform  input. 
This  corresponds  to  a  frequency  of  7.85  MHz.  Channel  2  has  no  direct  input  signal, 
but  it  measures  a  small  peak  to  peak  voltage  induced  from  the  adjacent  channel’s 
waveform. 
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4.5  Data  Collection  Components  and  Calibration 

For  experimental  operation,  the  PDA  collimator  and  amplifier  box  are  mounted  on 
the  ZaP  vacuum  tank’s  4-5/8  in  conflats  with  3  in  diameter  fused  silica  view-ports 
using  a  circular  aluminum  flange  with  a  wide  slit  opening  in  the  center,  as  shown 
in  Figure  4.16.  Section  4.6  includes  a  more  detailed  discussion  of  the  tank  windows 
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Figure  4.16:  Photodiode  arrays  mounted  on  the  ZaP  experiment.  Note  that  the 
vacuum  tank  extension  is  not  shown  on  the  schematic.  All  experimentation  was 
conducted  at  the  z  =  0  plane. 


and  inner  electrode  openings.  From  the  DB-25  connectors  on  the  amplifier  box, 
each  output  signal  travels  to  analog  digitizers  along  50  ft  of  7  strand,  20  gauge,  32 
twisted  pair  wire  with  both  individual  and  overall  shielding  and  grounding.  To  avoid 
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generating  ground  loops,  the  cable  is  grounded  on  one  side  to  the  DB-25  connectors 
on  the  PDA  box  and  again  at  the  entrance  of  the  Faraday  cage,  or  screen  room 
containing  the  digitizers.  The  signal  wire  is  converted  from  the  twisted  pair  wire  to 
Ethernet  cable  using  RJ-45  jacks.  These  are  then  converted  to  BNC  plugs  mounted 
on  an  aluminum  box  at  the  top  of  the  screen  room,  as  shown  in  Figure  4.17.  BNC 
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Figure  4.17:  Screen  room  box  connections.  The  box  contains  the  grounding  wire  to 
the  screen  room,  conversion  from  the  32  twisted  pair  cable  to  RJ-45  jacks  and  ethernet 
to  BNC  plugs. 
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to  LEMO  connector  cables  connect  this  box  to  two,  16  channel  Joerger  Model  TR 
analog  digitizers.  Each  12  bit  digitizer  operates  at  40  MHz  and  has  a  ±2  V  input 
range.  Unfortunately,  the  voltage  magnification  by  the  amplifier  boards  requires  the 
use  of  passive  signal  attenuation  to  reduce  the  digitizer  input  signal  to  the  range 
indicated.  10  and  14  dB  attenuators  placed  in-line  between  the  BNC  pins  on  the 
screen  room  box  and  the  BNC  to  LEMO  cables  as  shown  in  Figure  4.17  prevent  the 
input  signals  from  saturating  the  digitizers. 

Attenuation  factors  are  generated  by  measuring  a  known  waveform  with  a  1- 
2  MHz  frequency  through  each  attenuator  using  the  circuit  diagram  in  Figure  4.18. 
The  attenuation  value  is  simply  the  voltage  out  of  the  attenuator  divided  by  the 


Figure  4.18:  Attenuator  calibration  circuit  schematic.  In  figure  (a),  a  known  voltage 
Vo  is  applied  in  parallel  through  an  attenuator  and  a  1  MQ  resistor.  Figure  (b)  is  an 
enlarged  view  of  the  attenuator. 


voltage  in,  given  by  the  equation 
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where  resistor  R±  =  R2  and  follow  the  relationship 
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These  values  are  included  in  the  calibration  factors  used  for  data  analysis. 


(4.1) 


(4.2) 
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4.6  Viewing  Constraints  and  Window  Calibration 

Variations  in  the  vacuum  tank  window  and  inner  electrode  opening  size  also  affect  the 
signal  detected  by  the  photodiodes.  As  shown  in  Figure  4.19,  the  two  PDAs  view  the 
plasma  through  a  3.0  in  diameter  fused  silica  window  and  a  1.86  in  diameter  outer 
electrode  hole.  This  small  opening  restricts  the  viewing  axis  of  the  array  such  that 
the  1st  and  31st  collimator  holes  clip  the  sides  of  the  outer  electrode,  and  channel  32 
does  not  directly  view  the  opening.  Therefore,  the  data  taken  by  the  arrays  represent 
the  complete  observable  area  through  the  view-ports. 

The  fused  silica  vacuum  tank  window  also  proved  to  be  a  significant  contributor 
to  variations  in  the  collected  data.  During  the  course  of  experimentation,  the  plasma 
deposits  tungsten  metal  and  other  various  elements  onto  the  view-ports.  The  accrue¬ 
ment  of  this  material  distorts  the  signal  seen  by  the  PDAs  by  diffusing  the  plasma 
emission  through  the  window.  Attempts  to  calibrate  these  distortions  proved  incon¬ 
clusive  and  led  to  the  removal  and  sonic  cleaning  of  the  fused  silica.  The  Strobotac® 
strobe  light  used  for  channel  calibration  was  mounted  adjacent  to  the  view-port  as 
shown  in  Figure  4.20  and  set  to  high  luminosity  with  a  flash  frequency  of  5  Hz.  A 
frosted  glass  window  was  also  placed  between  the  strobe  and  the  view-port  opening 
to  achieve  a  more  uniform  intensity,  and  the  Joerger  digitizers  collected  data  from  a 
single  light  pulse  by  using  an  optical  cable  to  trigger  a  measurement.  Averaged  plots 
of  these  data  through  the  old  and  new  window,  as  well  as  through  no  window  at  all 
are  presented  in  Figure  4.21. 

The  plots  in  this  figure  show  significant  diffusion  of  the  strobe  light  through  the 
old  windows  that  would  make  any  sort  of  relative  comparison  between  the  two  PDA’s 
questionable.  While  the  new  window  and  no  window  data  are  similar  from  channel  to 
channel,  remaining  variations  in  the  channel  calibration  factors  and  frosted  glass,  non¬ 
uniformity  of  the  source  light,  and  misalignment  of  the  strobe  light  and  viewing  PDA 
preclude  this  relationship  from  being  used  to  generate  window  calibration  factors. 


Figure  4.19:  Viewing  area  of  the  Z-pinch  seen  by  two  PDAs  at  the  z  —  0  windows. 
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Figure  4.20:  Window  calibration  of  the  PDA  on  the  ZaP  tank.  The  photodiode  array 
being  tested  is  mounted  on  the  opposite  window  viewing  the  strobe  light.  The  fused 
silica  window  adjacent  to  the  strobe  is  removed  for  the  test,  and  data  are  taken  for 
the  dirty  window,  new  window,  and  no  window  on  the  PDA  side. 


Further  attempts  to  use  an  opal  glass  lens  to  better  diffuse  the  strobe  light  also 
proved  inconclusive.  In  the  end,  shots  taken  with  the  clean  windows  are  assumed  to 
have  no  attenuation  across  all  of  the  channels. 
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(a)  PDA  1  Window  Calibration 


-*=  PDA1  Old  Window  -"-PDAI  No  Window  PDA1  New  Window 


(b)  PDA  2  Window  Calibration 


Figure  4.21:  Averaged  frosted  glass  window  calibration  data  for  both  PDAs.  The 
frosted  glass  placed  between  the  strobe  and  the  tank  as  shown  in  Figure  4.20  produced 
a  centered  profile  instead  of  the  desired  uniform  profile.  Note  that  the  debris  build¬ 
up  on  the  old  windows  has  significantly  reduced  and  diffused  the  signal.  The  spike 
at  channel  22  was  the  result  of  an  improperly  terminated  digitizer  channel.  Other 
channel  variations  are  primarily  due  to  amplifier  calibration  problems  and  variations 
in  the  frosted  glass.  Without  a  true  uniform  profile,  it  is  not  possible  to  produce 
accurate  window  calibration  factors. 
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Chapter  5 

EXPERIMENTAL  RESULTS 

After  cleaning  the  vacuum  tank  windows,  all  emissivity  data  from  the  PDAs  were 
collected  at  the  z  —  0  plane  indicated  in  Figure  4.16.  A  bank  charge  of  5-9  kV  was 
utilized  for  all  shown  data.  The  majority  of  the  data  listed  in  this  chapter  are  from 
the  5  kV  shots,  where  more  useful  data  was  available  due  to  a  longer  quiescent  period. 
Ha  and  bremsstrahlung  filters  were  also  used  to  view  different  regions  of  the  spectra. 
The  bremsstrahlung  filter  in  particular  allows  for  the  comparison  of  PDA  emissivity 
against  chord-integrated  electron  density  measured  by  the  interferometer,  as  previ¬ 
ously  discussed  in  Section  3.4.  Comparisons  between  measured  Imacon  and  PDA 
emissivity,  and  magnetic  center  of  current  and  emission  centroid  are  also  determined 
from  these  chord-integrated  data. 

An  axisymmetric,  piecewise  constant  radial  profile  of  the  ZaP  Flow  Z-Pinch  was 
determined  from  the  1-D  Abel  inversion  technique  described  in  Section  3.6.  The  2-D 
Fourier  deconvolution  technique  from  Section  3.9  provides  a  more  detailed  view  of  the 
plasma,  with  both  a  radial  and  azimuthal  piecewise  constant  profile.  This  technique 
also  tracks  the  movement  of  the  profile  off  the  z-axis  and  can  be  used  to  investigate 
variations  between  the  emissions  of  unfiltered,  Ha ,  and  bremsstrahlung  filtered  light. 
Discussions  of  these  results  are  included  in  this  chapter. 

5.1  Chord-Integrated  Results 

PDA  emissivity  data  stored  on  the  Joerger  Model  TR  digitizers  were  run  through  the 
MATLAB  codes  using  the  theory  described  in  Chapter  3.  Each  PDA  used  16  channels 
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to  collect  data  through  the  x-  and  y- axes  of  the  tank  at  the  z  =  0  plane.  The  ma¬ 
jority  of  the  channels  are  evenly  spaced;  however,  problems  resulted  when  variations 
between  channel  calibration  factors  resulted  in  signals  with  disproportionably  larger 
or  smaller  gain  than  their  neighbors.  Studying  the  raw  output  data  revealed  where 
these  problems  were  pronounced,  and  care  was  taken  to  take  emission  measurements 
with  channels  where  these  signal  variations  were  minimal.  Unfortunately,  a  bad  at¬ 
tenuator  on  the  screen-room  box  shown  in  Figure  4.17  resulted  in  the  removal  of  one 
channel  from  the  x-axis  viewing  PDA. 

Magnetic  data  taken  for  a  5  kV  shot  shows  a  defined  current  sheet  during  plasma 
formation  and  an  approximately  1  cm  radius  peaked  emission  profile  during  the  quies¬ 
cent  period  described  by  the  magnetic  probes.  These  results  are  shown  in  Figure  5.1. 
The  unfiltered  emission  also  shows  movement  of  the  pinch  along  each  independent 
viewing  axis,  as  well  as  pronounced  amplitude  “flashes”  during  quiescence. 

Filtered  data  also  provide  useful  information  about  the  ZaP  Flow  Z-Pinch.  Filter¬ 
ing  the  viewed  spectra  involved  simple  Kodak-Wratten  # 12  and  #58  filters.  These 
him  filters  are  placed  between  the  collimator  top  plate  and  the  vacuum  tank  window. 
The  line  radiation  detected  by  the  photodiodes  previously  shown  in  Figure  4.3(b) 
shows  that  N-III  is  the  primary  contributor  to  emissivity.  These  lines  are  further 
reduced  by  the  filters  according  to  the  response  curves  shown  in  Figure  5.2.  The  #12 
filter,  (a)  and  (b),  only  allows  transmission  of  radiation  above  510  nm,  effectively 
eliminating  all  major  lines  other  than  Ha.  Ha  is  of  interest  because  it  represents  the 
emission  given  off  due  to  collisions  between  neutral  gas  in  the  plasma  with  hotter  ions. 
This  sort  of  interaction  occurs  at  the  boundary  of  the  pinch,  therefore  a  emission  fil¬ 
tered  for  Ha  should  have  a  hollow  profile  at  the  pinch  radius.  The  #58  is  a  band-pass 
filter  that  only  passes  radiation  between  500-580  nm.  When  combined  with  the  #12 
filter,  only  radiation  between  510-580  nm  is  detected  by  the  photodiodes.  This  region 
of  the  ZaP  plasma  has  very  little  line  radiation,  therefore  combining  the  filters  allow 
the  PDAs  to  measure  primarily  continuum  bremsstrahlung  radiation. 
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Normalized  m  =  1  Magnetic  Mode  Displacement 


Figure  5.1:  Unfiltered  chord-integrated  results  for  2  PDAs  viewing  a  5  kV  bank  voltage 
plasma.  Peaked  profiles  with  an  approximately  1  cm  radius  are  observed  during  the 
quiescent  period  described  by  the  normalized  m  =  1  mode  data  and  shown  between 
the  vertical  dashed  lines  on  all  plots.  Pinch  movement  along  the  viewing  axis  is  also 
observed  during  this  time,  as  indicated  by  the  centroid  position  shown  in  red.  During 
the  formation  period,  the  PDAs  record  the  current  sheet  traversing  the  z  =  0  plane, 
shown  in  this  pinch  at  approximately  24-27  /zs. 
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ZaP  Plasma  Emission  Spectra 
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Figure  5.2:  Filter  curves  and  viewed  spectra  of  the  ZaP  Flow  Z-Pinch.  (a)  The  blue 
line  is  the  spectra  seen  by  the  photodiodes  without  filtering.  The  red  line  repre¬ 
sents  the  spectral  response  curve  of  the  Kodak-Wratten  # 12  filter,  (b)  This  filter 
leaves  Ha  as  the  primary  contributor  to  line  radiation.  (c,d)  The  response  curve  and 
viewed  spectra  for  the  ^58  filter  still  leaves  some  line  radiation  around  500  nm.  (e) 
Combining  the  two  filters  produces  a  narrower  response  curve  between  510-580  nm, 
thereby  removing  almost  all  line  radiation  from  the  ZaP  plasma,  (f)  This  leaves  only 
bremsstrahlung  radiation  as  the  primary  emission  seen  by  the  PDAs. 
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Comparisons  between  5  kV  bank  voltage  shots  for  these  filters  are  shown  in  Fig¬ 
ure  5.3.  Although  the  filters  significantly  reduce  the  brightness  of  the  measured  light, 
especially  in  the  bremsstrahlung  case,  the  PDA  is  still  able  to  detect  a  peaked  profile 
during  quiescence.  These  data  are  therefore  useful  for  further  analysis  later  in  this 
chapter. 

5.2  Electron  Density  Results  and  Comparisons 

Evaluation  of  the  emission  profile  from  bremsstrahlung  radiation  and  data  from  the 
interferometer  shows  that  the  electron  density  of  the  plasma  is  proportional  to  the 
square  root  of  the  emissivity.  Interferometry  data  taken  through  the  y-axis  at  z  —  0 
prior  to  mounting  the  side-mounted  PDA  provides  line-integrated  electron  density 
of  the  plasma  along  a  chord  through  the  center  of  the  window  at  z  =  0  and  along 
a  second  chord  1  cm  off  axis.  These  data  show  a  peaked  density  plasma  profile 
during  quiescence,  with  higher  density  measured  along  the  center  chord.  Although 
the  interferometer  is  removed  to  allow  both  PDAs  to  view  the  plasma  at  z  =  0, 
relating  5  kV  bank  voltage  shots  with  similar  m  —  1  magnetic  mode  data  allows 
for  general  comparison  between  the  interferometry  and  PDA  signals.  As  described 
in  Section  3.4,  bremsstrahlung  emissivity  E  oc  n2e.  Figure  5.4  shows  that  the  signal 
detected  by  PDA  channels  at  similar  centered  and  offset  locations  with  respect  to  the 
x-  and  y- axes  also  display  this  same  peaked  profile.  Since  the  PDA  and  interferometer 
are  looking  at  different  plasmas,  amplitude  variations  in  the  two  chords  indicative  of 
pinch  movement  are  not  similar;  however,  the  signal  at  the  center  chord  location 
generally  remains  higher  than  the  offset  channel  during  quiescence. 

5.3  Visual  Emission  Results  and  Comparisons 

The  Imacon  streak  photos  previously  discussed  in  Section  2.4  provide  an  independent 
verification  of  the  photodiode  emissivity.  Since  the  time  scale  for  the  Imacon  photo  is 
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Figure  5.3:  Chord-integrated  profiles  of  (a)  unfiltered,  (b)  Ha,  and  (c)  bremsstrahlung 
filtered  5  kV  bank  voltage  plasma  as  seen  by  a  16  channel  PDA.  The  vertical  yellow 
lines  represent  the  boundaries  of  the  quiescent  period.  The  red  line  is  the  calculated 
emissivity  centroid  of  the  plasma.  Note  that  the  measured  signal  drops  with  increased 
filtering;  however,  a  peaked  profile  is  still  evident  in  all  cases  during  quiescence. 
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Figure  5.4:  Bremsstrahlung  emission  and  line-integrated  electron  density,  (a)  shows 
the  line-integrated  electron  number  density  measured  by  the  interferometer  for  a  5  kV 
bank  voltage  plasma  at  a  center  (black  line)  and  1  cm  offset  (red  line),  (b)  and  (c) 
are  the  square  roots  of  the  chord-integrated  emission  profiles  for  a  bremsstrahlung 
filtered  5  kV  bank  voltage  plasma  at  centered  and  approximately  1  cm  offset  channels 
along  the  vertical  and  horizontal  axis  respectively.  The  legends  indicate  the  positions 
of  the  chords.  Although  the  measurements  are  for  different  pinches,  the  shots  exhibit 
similar  (d)  magnetic  mode  data  and  quiescent  period  boundaries,  shown  as  the  blue 
and  green  dashed  lines.  Both  the  electron  number  density  and  the  square  root  of 
emission  show  a  peaked  plasma  profile  on  axis  during  quiescence  with  a  lower  density 
or  emissivity  at  a  1  cm  radius. 


70 


roughly  500  ns  per  millimeter  and  runs  approximately  8  cm  in  length,  beginning  the 
image  during  the  formation  stage  allows  for  the  majority  of  the  quiescent  period  to  be 
recorded.  Mounted  opposite  the  y-axis  viewing  PDA  at  z  —  0,  the  chord-integrated 
emissivity  and  vertical  movement  recorded  by  both  diagnostics  during  quiescence 
is  very  similar.  The  centroid  of  emissivity  determined  using  the  theory  discussed 
in  Section  3.2  also  correlates  to  the  camera  emission  recorded  on  the  film.  Slight 
differences  in  signal  amplitude  are  most  likely  the  result  of  image  saturation  and  film 
sensitivity  of  the  Irnacon  photos.  Variations  in  the  amplifier  calibrations  and  scaling 
factors  also  result  in  differences  between  the  profiles.  Figure  5.5  shows  similarities 
between  the  PDA  centroid  and  Imacon  emission  profile  for  an  8  kV  bank  voltage 
plasma  during  a  short  quiescent  period  between  29-36  [is. 

Also  shown  in  Figure  5.6  is  a  5  kV  bank  voltage  shot  with  a  much  longer  qui¬ 
escent  period.  Both  of  the  photos  and  chord-integrated  data  are  unfiltered,  and  the 
Imacon  image  is  fit  to  a  y-axis  similar  to  the  one  set  for  the  PDA  profile.  With  these 
estimations,  the  PDA  and  Imacon  contours  show  the  same  y-axis  movement  during 
formation  and  amplitude  variations  during  quiescence. 


5.4  Centroid  Results  and  Comparisons 

Analysis  of  the  center  of  current  measurements  using  the  magnetic  probes  described  in 
Section  2.2  and  the  emissivity  centroid  discussed  in  Section  3.2  confirms  that  a  direct 
relationship  between  the  two  exists  during  the  majority  of  quiescence.  Figure  5.7 
shows  an  unfiltered  5  kV  bank  voltage  plasma  viewed  in  the  x-axis  by  15  channels 
of  a  PDA  and  an  8  probe  azimuthal  magnetic  probe  array  at  the  z  —  0  plane.  The 
azimuthal  magnetic  probes  measure  both  the  amplitude  and  phase  angle  associated 
with  the  rn  —  1  mode.  Projecting  these  data  onto  the  x-axis  and  multiplying  it 
by  the  wall  radius  divided  by  two  gives  the  x-axis  center  of  current.  In  general, 
both  diagnostic  systems  assign  similar  locations  for  emission  and  current  centers. 
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Figure  5.5:  (a)  Imacon  and  (b)  PDA  views  of  an  unfiltered  8  kV  bank  voltage  plasma. 
Both  diagnostics  view  the  pinch  through  the  z  —  0  plane’s  side  windows.  Note 
that  both  systems  record  the  vertical  pinch  movement.  The  red  line  on  both  plots 
represents  the  PDA  calculated  centroid  position. 
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Figure  5.6:  (a)  Imacon  and  (b)  PDA  views  of  an  unfiltered  5  kV  bank  voltage  plasma 
showing  pinch  quiescence.  Imacon  film  saturation  and  PDA  emission  scaling  result 
in  the  different  emissivity  amplitudes.  The  red  line  in  both  plots  represents  the  PDA 
calculated  centroid  position. 
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Figure  5.7:  Emission  centroid  and  magnetic  center  of  current  for  an  unfiltered  5  kV 
bank  voltage  Z-pinch.  The  red  line  represents  the  PDA  calculated  emissivity  centroid. 
The  yellow  line  is  the  x-axis  component  of  the  rri  —  1  mode  data  representing  the 
center  of  current  measured  by  the  azimuthal  magnetic  probes.  The  vertical  yellow 
lines  delineate  the  quiescent  period  boundaries. 
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Variations  between  the  measurements  are  due  problems  with  the  magnetic  probe  array 
as  well  as  minor  calibration  inaccuracies  in  the  PDA  amplifiers.  As  the  magnetic  probe 
data  cannot  fully  resolve  higher  frequency  bulk  plasma  motion,  any  quick  movement  of 
the  emission  center  is  not  accurately  followed  by  the  magnetic  probe  data;  however, 
general  comparison  shows  a  direct  correlation  between  the  emissivity  centroid  and 


current  center. 
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5.5  Abel  Inversion  Results 

Axisymmetric  1-D  profiles  determined  using  an  Abel  inversion  show  flat  radial  pro¬ 
files  during  pinch  formation,  and  a  peaked  radial  pinch  with  a  1  cm  radius  during 
quiescence.  The  increasing  disparities  between  plasma  emission  center  and  centroid 
position  also  play  an  important  role  in  profile  construction.  As  the  plasma  emission 
center  and  centroid  diverge,  differences  between  the  left  and  right  half  path-length 
matrices  increase  in  magnitude,  causing  discontinuities  between  the  center  locations 
of  each  side.  Figure  5.8  is  an  Abel  profile  of  unfiltered  5  kV  bank  voltage  plasma 
determined  by  a  15  chord  of  a  PDA  viewing  along  the  y-axis  of  the  pinch.  During 
formation,  the  normalized  chord-integrated  profile  remains  peaked  while  the  radial 
profile  remains  flat.  In  quiescence,  the  radial  profile  exhibits  a  peaked  contour  with 
an  approximately  1  cm  radius.  A  discontinuity  arises  between  the  centers  of  each 
side  of  the  radial  emissivities  when  the  centroid  algorithm  cannot  accurately  follow 
the  plasma  contour  in  an  off-center  profile.  However,  when  the  emission  centroid  ac¬ 
curately  follows  the  plasma  center,  the  inverted  profile  exhibits  continuity  across  the 
left  and  right  sides.  In  general,  since  the  pinch  remains  near  the  center  during  quies¬ 
cence,  the  inversion  profile  at  this  time  accurately  models  the  peaked  chord-integrated 
contour. 

5.6  Fourier  Deconvolution  Results 

The  2-D  Fourier  technique  described  in  Section  3.10  provides  a  piecewise  constant  ra¬ 
dial  and  cosine  varying  azimuthal  plasma  profile  of  the  ZaP  Z-pinch  using  the  theory 
from  Section  3.10.  Data  from  unfiltered,  Ha  and  bremsstrahlung  filtered  emissions 
all  show  movement  of  the  pinch  profile  off  axis  and  amplitude  variations  during  qui¬ 
escence.  Two  16  channel  PDAs  mounted  at  z  —  0  collected  chord-integrated  data  as 
previously  mentioned  in  Section  5.1.  The  MATLAB  code  Full2d.m  takes  these  data 
and  produces  a  two-dimensional  plasma  centered  array  of  emissivity  values  versus 
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Figure  5.8:  Abel  inverted  axisymmetric  radial  profiles  for  an  unfiltered  5  kV  bank  volt¬ 
age  Z-pinch.  (a-d)  The  black  contour  line  represents  the  normalized  chord-integrated 
signals  along  the  y-axis.  The  red  and  blue  lines  are  the  inverted  radial  profiles  for  the 
left  and  right  half  sides.  The  points  represent  the  discrete  data  points.  The  dashed 
yellow  lines  in  the  y-axis  contour  plot  show  the  boundaries  of  the  quiescent  period 
while  the  solid  cyan  lines  show  the  time  of  each  shown  inverted  profile,  (a)  During 
formation,  the  radial  profile  of  the  pinch  remains  flat,  (b)  In  quiescence  however, 
it  has  a  peaked  radial  profile  with  a  radius  of  about  1  cm.  (c)  When  the  centroid 
and  emission  profile  do  not  accurately  match;  discontinuity  arises  between  the  left 
and  right  sides  of  the  radial  profile,  (d)  In  cases  where  the  centroid  and  emission  do 
match  however,  the  radial  profile  remains  continuous  across  the  center. 
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time.  A  separate  vector  of  x  and  y  centers  chosen  by  the  TDerror  .m  and  bestf  it  .m 
functions  combines  with  the  array  data  to  provide  lab-frame  profiles  of  the  Z-pinch 
versus  time. 

Figure  5.9  shows  the  x-  and  y-axis  chord-integrated  profiles  of  an  unfiltered  5  kV 
bank  voltage  plasma  while  Figure  5.10  shows  several  deconvolved  profiles  of  this  pinch 
at  various  times  during  quiescence.  Peaked  and  hollow  profiles  of  the  pinch  have  a 
radius  of  about  1  cm  and  exhibit  varying  position  and  amplitude.  During  quiescence, 
movement  off  axis  is  minimal  and  changes  in  emission  amplitude  are  smooth.  In 
contrast,  during  the  post  quiescent  period  the  plasma  movement  and  amplitude  both 
change  wildly. 

Figures  5.11  and  5.12  are  of  a  similar  5  kV  bank  voltage  plasma  viewed  through 
the  Kodak- Wratten  #12  filter  to  view  Ha  emissions.  This  plasma  exhibits  a  hollow 
profile  throughout  the  quiescent  period  with  peaks  along  a  radius  of  about  1  cm  from 
the  center  of  the  profile.  This  is  indicative  of  the  Ha  shell,  where  cold  neutral  gas  hit 
by  the  plasma  emits  the  majority  of  the  radiation  seen  through  the  filtered  lens.  As 
expected,  the  amplitude  of  this  emission  is  lower  than  the  unfiltered  case. 

Figures  5.13  and  5.14  are  the  chord-integrated  profiles  and  2-D  Fourier  profiles 
of  a  5  kV  bank  voltage  plasma  viewed  through  the  bremsstrahlung  filters  listed  in 
Section  5.1.  As  in  the  unfiltered  case,  both  peaked  and  hollow  profiles  are  seen  in 
the  deconvolved  radial  contours,  however  the  signal  amplitudes  are  much  lower.  The 
radius  of  these  profiles  is  also  about  1  cm  in  radius  and  moves  off  the  2-axis  during 
quiescence.  During  periods  of  low  chord-integrated  signal  in  the  quiescent  period,  the 
2-D  cmissivity  profile  exhibits  a  hollow  shape,  as  seen  in  Figure  5.14(b).  When  the 
signal  increases  during  a  “flash,”  emission  in  the  center  of  the  2-D  contour  increases, 
resulting  in  a  peaked  profile  as  shown  in  Figures  5.14(a),(b),(d). 
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Figure  5.9:  x-  and  y-axis  contours  of  an  unfiltered  5  kV  bank  voltage  Z-pinch.  16 
chords  view  the  plasma  along  the  y-axis  and  15  view  it  along  the  x-axis.  The  dashed 
yellow  lines  represent  the  boundaries  of  the  quiescent  period.  The  solid  blue  lines  are 
the  times  the  2-D  data  shown  in  Figure  5.10  were  taken. 
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Figure  5.10:  2-D  emissivity  profiles  of  the  unfiltered  5  kV  bank  voltage  shown  in 
Figure  5.9.  (a)  and  (b)  show  a  partially  hollow  profile.  The  axisymmetric  inner  shell 
is  a  product  of  the  constraint  put  on  the  Fourier  deconvolution  analysis  discussed  in 
Section  3.11.  (c)  shows  a  peaked  pinch  profile  with  a  1  cm  radius,  (d)  is  a  peaked 
profile  with  a  1  cm  offset  in  the  positive  ic-axis.  This  occurs  immediately  before  the 
pinch  goes  unstable,  as  shown  by  the  last  vertical  blue  line  in  Figure  5.9.  Also  note 
that  the  shells  chosen  for  this  deconvolved  profile  have  larger  radii  and  coarser  spacing 
due  to  the  offset  reducing  the  number  of  chords  viewing  the  center  of  the  plasma. 
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Figure  5.11:  Ha  chord-integrated  x-  and  y-axis  contours  of  a  5  kV  bank  voltage  Z- 
pincli.  16  chords  view  the  plasma  along  the  y- axis  and  14  view  it  along  the  x-axis. 
The  dashed  yellow  lines  represent  the  boundaries  of  the  quiescent  period.  The  solid 
blue  lines  are  the  times  2-D  data  shown  in  Figure  5.12  were  taken. 
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Figure  5.12:  2-D  profiles  of  Ha  emission  from  the  5  kV  profile  in  Figure  5.11.  (a-d) 
Throughout  the  quiescent  period  of  the  pinch,  the  2-D  radial  profile  remains  hollow 
with  a  radius  of  approximately  1  cm.  This  is  due  to  collisions  between  the  plasma 
and  colder  neutral  gas  at  that  radius  radiating  the  majority  of  the  Ha  radiation. 
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Figure  5.13:  Bremsstrahlung  chord-integrated  x-  and  y-axis  contours  of  a  5  kV  bank 
voltage  Z-pinch.  16  chords  view  the  plasma  along  the  y-axis  and  15  view  it  along 
the  x-axis.  The  dashed  yellow  lines  represent  the  boundaries  of  the  quiescent  period. 
The  solid  blue  lines  are  the  times  2-D  data  were  taken  in  Figure  5.14. 
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Figure  5.14:  2-D  profiles  of  bremsstrahlung  emission  from  the  5  kV  profile  in  Figure 
5.13.  (a),  (c),  and  (d)  show  a  peaked  profile  with  a  1  cm  radius  during  periods  of 
bright  emission  flashes,  (b)  When  the  emission  is  lower  during  the  quiescent  period, 
such  as  indicated  by  the  third  blue  line  in  Figure  5.13,  the  2-D  profile  has  a  hollow 
shape.  Overall  emission  is  also  lower  compared  to  the  un filtered  and  Ha  case. 
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Chapter  6 

CONCLUSIONS 

Emissivity  profiles  and  tomographic  data  taken  by  both  PDAs  prove  valuable  in 
understanding  the  ZaP  Flow  Z-Pinch.  Using  the  theory  and  diagnostics  discussed  in 
previous  chapters  allows  for  the  comparison  of  PDA  emissivity  data  to  measurements 
made  by  the  Irnacon  camera,  magnetic  probes,  and  interferometer.  Verification  of 
other  diagnostic  systems  by  the  PDAs  increases  the  confidence  in  data  interpretation. 
Radial  1-D  Abel  inverted  and  2-D  Fourier  deconvolved  profiles  of  emissivity  are  also 
constructed  from  PDA  data,  giving  a  better  understanding  of  the  emissivity  profile 
of  the  ZaP  Flow  Z-Pinch.  Several  simple  improvements  and  alternative  experimental 
setups  are  also  presented  in  this  chapter. 

6.1  Summary 

Experimental  results  of  chord  integrated  data  show  that  the  emissivity  profile  detected 
by  the  PDAs  is  peaked  during  quiescence  and  matches  the  visual  profile  measured 
by  the  Imacon  streak  photos.  The  PDAs  accurately  measure  amplitude  variations 
and  movement  of  the  pinch  along  the  viewing  axis.  Ha  and  bremsstrahlung  filters 
significantly  decrease  the  detected  signal;  however,  the  data  still  show  features  of  these 
types  of  plasma  emissions  during  quiescence  and  are  therefore  useful  for  analysis. 

Light  filtered  for  bremsstrahlung  radiation  for  5  kV  bank  voltage  pinches  show 
a  peaked  profile  on  axis  similar  to  interferometry  measurements  of  line  integrated 
electron  density  of  5  kV  shots  with  comparable  normalized  mode  data.  This  provides 
a  general  validation  that  bremsstrahlung  radiation  is  proportional  to  the  square  root 
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of  the  electron  number  density. 

The  centroid  of  emissivity  determined  from  the  centroid. m algorithm  also  relates 
to  the  center  of  current  determined  from  the  azimuthal  magnetic  probes.  In  general, 
the  emission  and  the  centroid  appear  to  overlay  during  quiescence.  As  the  emission 
moves  further  off  the  2- axis  however,  the  centroid  data  begins  to  diverge  from  the 
viewed  peaked  emission.  This  validates  the  truncation  of  these  points  to  reduce 
errors  in  the  Abel  and  Fourier  algorithms. 

The  Abel  inversion  algorithm  Abel  .m  correctly  deconvolves  properly  centered  ax- 
isymmetric  synthetic  and  experimental  radial  profiles.  Results  from  5  kV  bank  voltage 
plasma  show  that  the  Z-pinch  maintains  a  1  cm  peaked  and  symmetric  profile  during 
quiescence.  When  the  center  of  the  inversion  profile  significantly  differs  from  the  cen¬ 
ter  of  the  chord  integrated  contour,  discontinuities  arise  between  the  center  channels 
of  the  left  and  right  half  sides. 

2-D  Fourier  deconvolved  data  determined  using  the  Full2d.m program  accurately 
measures  experimental  and  synthetic  movement  off  the  z-axis,  amplitude  variations 
in  emission,  and  both  hollow  and  peaked  plasma  contours.  Experimental  5  kV  bank 
voltage  plasmas  exhibit  a  radius  of  about  1  cm  during  quiescence  regardless  if  they 
are  hollow  or  peaked.  The  TDerror.m  and  bestfit.m  algorithms  accurately  offset 
the  shell  center  and  change  shell  radii  to  produce  the  best  least  squares  fit  of  the 
given  data.  Unfiltered  data  have  the  largest  amplitude  and  a  combination  of  hollow 
and  peaked  profiles  during  quiescence.  Ha  emissions  viewed  through  a  filter  show  a 
hollow  profile  with  a  radius  of  1  cm,  corresponding  to  the  estimated  Ha  shell  size  of 
the  Z-pinch.  Light  filtered  for  bremsstrahlung  shows  hollow  and  peaked  profiles  with 
a  lower  signal.  During  times  of  greater  intensity,  the  emissivity  from  this  radiation 
appears  to  be  more  peaked  than  at  other  times  in  the  quiescent  period,  where  it 
exhibits  a  more  hollow  shape. 


85 


6.2  Future  Work 

Although  the  photodiode  arrays  already  provide  a  great  deal  of  information  regarding 
the  ZaP  Flow  Z-Pinch  plasma,  many  simple  calibrations  and  changes  to  experimental 
setup  can  improve  the  quantity  and  quality  of  the  data  collected.  Problems  with 
amplifier  calibrations,  vacuum  tank  windows,  cabling,  and  attenuators  continuously 
plagued  data  collection  and  skewed  data  interpretation.  Although  the  current  design 
is  much  more  robust  than  the  earlier  prototype  PDA,  the  electronic  circuits  were  not 
completely  reliable. 

Vibrations  from  mounting  and  handling  the  PDA,  as  well  as  those  from  the  ex¬ 
periment  itself  caused  circuit-board  connections  to  come  loose,  causing  changes  in  the 
amplifier  responses.  Sealing  the  trimmer  capacitors  on  the  amplifier  boards  shown  in 
Figures  4.8  and  4.10  with  an  epoxy  resin  should  prohibit  further  damping  changes  in 
the  PDA  response.  Once  completed,  new  calibration  factors  should  be  determined 
using  the  process  described  in  Section  4.3.  The  weight  of  the  36  twisted  pair  shielded 
cable  also  puts  a  great  deal  of  strain  on  the  DB-25  connectors  shown  in  Figure  4.11. 
Proper  strain  relief  has  reduced  this  force,  however  resoldering  of  the  connector  pins 
is  required  to  reduce  loose  channel  connections  that  appear  occasionally  during  ex¬ 
perimentation. 

Better  calibration  of  the  vacuum  tank  windows  will  give  a  much  more  accurate 
view  of  the  plasma  profile.  A  constant  light  source  with  a  flat  profile  over  the  entire 
viewing  area  is  required  to  accurately  measure  window  variations  and  dispersion. 
Combining  several  opal  glass  plates  will  properly  flatten  a  light  source  but  will  also 
require  a  high  intensity  light  to  make  up  for  the  reduction  in  signal  through  each 
plate. 

Reducing  variations  between  the  center  of  current  and  centroid  of  cmissivity  will 
require  new  calibration  factors  for  the  magnetic  probes.  This  recalibration  is  currently 
underway  on  the  experiment,  and  simply  replacing  the  calibration  values  should  re- 
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move  many  of  the  discrepancies.  In  addition,  the  centroid  finder  routine  can  be 
improved  to  better  handle  Z-pinch  movement  further  off  the  center  axis. 

Although  this  thesis  presented  comparisons  between  line-integrated  electron  den¬ 
sity  measured  by  the  interferometer  and  bremsstrahlung  emission  viewed  by  the  PDA 
for  similar  plasmas  in  Section  5.1,  no  data  currently  exist  where  both  diagnostics  view 
the  same  local  area  of  the  plasma.  Removing  the  side  mounted  PDA  at  the  z  —  0 
plane  to  permit  the  interferometer  lasers  to  measure  y-axis  density  will  allow  for  this 
condition.  Although  the  remaining  PDA  measures  emission  along  the  re-axis,  this 
should  still  give  a  general  relationship  between  bremsstrahlung  emission  and  density 
for  the  same  pulse. 

Improvements  in  the  2-D  algorithm  performance  will  also  decrease  inaccuracies 
in  the  resulting  deconvolved  data.  The  2-D  profile  center  finder  TDerror .  m  does  not 
take  the  uncertainty  of  the  y2  position  into  account.  While  this  method,  discussed  in 
Section  3.11,  is  satisfactory  for  this  experimentation,  a  more  accurate  profile  centering 
can  be  attained  by  fitting  the  profile  according  to  the  variance  given  in  the  following 
equation, 

[Set]  =  [M2]  [5E2]  +  [8M2]  [E2]  .  (6.1) 

5e2  is  the  variance,  or  square  of  the  standard  deviation  of  the  radial  profile.  M  is  the 
position  of  y2  while  SM 2  is  the  uncertainty  of  that  measurement.  E  and  8E2  are  the 
chord  integrated  signals  and  signal  uncertainty  respectively.  The  current  algorithm 
only  checks  the  variance  associated  with  the  first  term  on  the  right  side.  By  including 
the  second  term,  a  more  accurate  fit  can  be  expected. 

Up  to  this  point,  all  data  taken  with  the  PDAs  have  been  in  one  plane  of  the 
z-axis.  Rotating  the  PDAs  so  that  the  chords  are  oriented  along  the  z-axis  will  allow 
for  comparison  between  chord-integrated  emissivity  and  the  wavelength  of  m  =  0  and 
m  —  1  perturbations  during  quiescence.  Listings  of  the  chord-integrated  codes  used 
in  data  analysis  of  the  PDA  are  given  in  Appendix  A. 
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Appendix  A 

MATLAB  CHORD  INTEGRATING  AND  ABEL 
INVERSION  CODES 

Appendix  A  includes  the  MATLAB  code  and  calling  functions  dealing  with  chord- 
integrated  and  Abel  inversion  data.  Explanations  of  each  code  are  listed  in  their  ap¬ 
propriate  subsection.  Codes  with  movie  or  plotting  routines  are  not  included  in  these 
listings.  Shown  below  is  the  function  calling  sequence.  The  overall  code  PDAcode.m 
takes  data  from  the  MDSplus  storage  system  and  creates  array  variables  of  the  chord 
integrated  signals  and  their  time  bases.  All  the  subsequent  routines  handle  data  from 
both  arrays  simultaneously,  with  the  exception  of  the  Abel  1 .  m  and  Abel2 .  m  functions. 


CALLING  SEQUENCE 

•  PDAcode.m 

—  tankaxis.m 

—  magneticmode. m 

—  centroid,  m 

—  centroiderror.m 

—  abell.m 


abel2.m 
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A.l  PDAcode.m 


7.===================================================================== 

°/o 

70Adam  Mads  on 

7.11/04/04 

% 

7.PH0T0DI0DE  ARRAY  EMISSIVITY  CALCULATIONS  PROGRAM:  PDAcode.m 

% 

7.Description:  This  is  the  overall  calling  program  used  for  analysis  of 
°/„  chord  integrated  and  Abel  inverted  PDA  data.  Information  is 
7.  called  from  the  MDS  plus  tree  using  the  mdsopen.m  and  mdsvalue.m 
7.  functions.  Commenting  out  various  call  functions  listed  below 
7.  will  allow  for  thecollection  of  the  data  of  interest  at  a  faster 

7.  rate.  The  dataplots  .mfunction  plots  out  all  chord  integrated, 

7.  centroid,  and  windowvalidation  plots  and  stores  them  as  .  eps  or 
%  -jPg  files. 

% 

/.===================================================================== 

°/o 

°/0begin  program 
clear  all;  close  all; 
format  short; 

7.===================================================================== 

7.7.7.  USER  INPUT  7.7.7. 

7.===================================================================== 

7.L00PING  SHOT  NUMBER 
7.shots  on  given  day 
7.  for  i  =  3:31 
7.  shotnum  =  50315000+i 

7.SINGLE  SHOT  NUMBER  CODES  (COMMENT  OUT  LOOP) 
shotnum  =  50316044 

7.Probe  Location 
probelocation  =  0; 

7.PR0BE  ORIENTATION 

7.  CHANGE  VALUES  IN  tankaxis.m  TO  REFLECT  PDA  ORIENTATION 


92 


’/.==================================================================== 

HI  MDS  DATA  ACCESS  V/X 

’/.==================================================================== 

“/POSSIBLE  LOOP  CREATION  TO  TRUNCATE  MDS  VALUE  BEFORE  PLACING  IN  ARRAY 

“/CODE  CHANGED  TO  ALLOW  FOR  PROCESSING  BOTH  PDA  DATA  SIMULTANEOUSLY 
“/„  -  PDA1  SPECIFIC  ARRAYS  WILL  END  IN  "1" 

l  -  PDA2  SPECIFIC  ARRAYS  WILL  END  IN  "2" 

°/„MDS  access  to  "alfven"  or  "landau"  computers 

“/„  [a,  status]  =  mdsopen( ’  landau .  aa.  Washington .  edu:  :  zapmain’  ,  shotnum) ; 
[a, status]  =  mdsopen( ’alfven. aa. Washington. edu: : zapmain’ , shotnum) ; 


’/.==================================================================== 

“/PDA  1  ARRAY  GENERATION 

PDA1  =  zeros  (32 , 16384)  ;  °/016384  from  16kB  digitzer  setting 

Tms_l  =  zeros (1 , 16384) ; 
nonzero  1  =  []  ; 

°/„Storing  MDS  data  to  Matlab  array 
for  i  =  1:32 

°/0PDAl  Frequency  Smoothed  Data 

channell  =  [’ of f set_baseline(smoothld’  ... 

’ (smoothld(\pdal_ ’ ,num2str (i)  ’,12), 13))’  ]; 

PDAl_tree  =  mdsvalue (channell) ; 

j  =  sum(size (PDAl_tree) ) ; 
if  j>1000 

PDAl(i,  :)  =  PDAl_tree’ ; 

end 

end 

j  =  0; 

for  i  =  1:32 

if  sum(abs (PDA1 (i , 1 : 500) ) )  ~=  0 

nonzero  1  (i-j )  =  i;  °/0index  of  row 

else 

j=j+i; 

end 


end 
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nonzero  1  °/„index  of  PDA  1  channels  with  collected  data 

"/PDA  use  check 

if  size (nonzerol)  ==  [0  0] 

PDAluse  =  0; 

else 

PDAluse  =  1; 

end 

if  PDAluse  ==  0 

XI  =  0;  time 1=0;  Tms_l=0;  PDA1=0;  ml=0;  nl=0; 
Max_Signall=0 ;  11=0;  Max_Signal_valuel=0; 

else 


XI  =  .  15748 .  *nonzerol ;  °/0 channel  position  (cm)  for  PDA1 

timel  =  [’dim_of  (\pdal_;>  ,num2str  (nonzerol  (1) )  ’)’]; 

Tms_l  =  mdsvalue (timel) 5 *le6 ;  "/time  in  microseconds 

"/SIGNAL  OUTPUT  FOR  PDA1 

PDA1  =  -PDA1  (nonzerol ,:)  ;  "/.returns  positive  signals 

"/Index  of  rows  and  columns 
7.  row:  channel,  column:  time  step 
[ml,nl]  =  size(PDAl); 

[Max_Signall ,  II]  =  max(PDAl);  °/Il  returns  max  signal  ch. 

Max_Signal_valuel  =  max(Max_Signall) ; 

end 


"/==================================================================== 

"/PDA  2  ARRAY  GENERATION 

PDA2  =  zeros  (32 , 16384) ;  °/016384  from  16kB  digitzer  setting 

Tms_2  =  zeros (1 , 16384) ; 
nonzero2  =  []  ; 

"/Storing  MDS  data  to  Matlab  array 
for  i  =  1:32 

°/0PDA2  Frequency  Smoothed  Data 

channel2  =  [)offset_baseline(smoothld)  ... 

’  (smoothld(\pda2_;’  ,num2str  (i)  ’,12),  13))"  ]; 
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PDA2_tree  =  mdsvalue (channel2) ; 

j  =  sum(size (PDA2_tree) )  ; 
if  j>1000 

PDA2 (i ,  :)  =  PDA2_tree’ ; 

end 

end 

j  =  0; 

for  i  =  1:32 

if  sum(abs (PDA2 (i , 1 : 500) ) )  ~=  0 

nonzero2(i-j )  =  i;  °/0index  of  row 

else 

j=j+i; 

end 

end 

nonzero2  7„index  of  PDA  2  channels  with  collected  data 

%PDA  use  check 

if  size (nonzero2)  ==  [0  0] 

PDA2use  =  0; 

else 

PDA2use  =  1; 

end 

if  PDA2use  ==  0 

X2  =  0;  time2=0;  Tms_2=0;  PDA2=0;  ra2=0;  n2=0; 

Max_Signal2=0 ;  12=0;  Max_Signal_value2=0 ; 

else 

X2  =  .  15748 .  *nonzero2 ;  °/0channel  position  (cm)  for  PDA2 

time2  =  [’dim_of (\pda2_’ ,num2str(nonzero2(l))  ’ )’]; 

Tms_2  =  mdsvalue (time2) J*le6;  %time  in  microseconds 

°/„SIGNAL  OUTPUT  FOR  PDA2 

PDA2  =  -PDA2 (nonzero2 , : ) ;  %returns  positive  signals 

°/0Index  of  rows  and  columns 
7.  row:  channel,  column:  time  step 
[m2,n2]  =  size(PDA2); 
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[Max_Signal2, 12]  =  max(PDA2);  °/„I2  returns  max  signal  ch. 

Max_Signal_value2  =  max(Max_Signal2) ; 

end 

“/WORKSPACE  CLEANING 

clear  timel  time2  status  i  j  channell  channel2  a  PDAl_tree  . . . 
PDA2_tree  II  12 

’/.===================================================================== 

HI  TANK  CENTERED  AXIS  HI 

’/.===================================================================== 

“/NOTE:  CHANGE  VARIABLES  » axisl ’  ,axis2)  TO  CORRESPOND  TO 
1  PDA  ORIENTATION 

[axisl , axis2]  =  tankaxis(nonzerol,nonzero2) ; 

“/simple  change  from  array  axis  to  tank  axis  channel  locations 
XI  =  axisl; 

X2  =  axis2; 

“/oWORKSPACE  CLEANING 
clear  axisl  axis2 


’/.==================================================================== 

HI  MAGNETIC  MODE  DATA  HI 

"/==================================================================== 

[magtime, mode, model, Cx,Cy, cutoff]  =  ... 
magneticmode (shotnum,probelocat ion) ; 

"/==================================================================== 

HI  ELECTRON  NUMBER  DENSITY  COMPARISON  HI 
’/.==================================================================== 

%Use  on  Brem  Filter  shots  and  compare  to  similar  pO  interferometer 
1  shots  E  proportional  ne~2 

1  interf erometer (PDA1 ,Tms_l ,X1 ,nonzerol , shotnum, mode, magtime, cutoff ) ; 


"/================================= 

HI  CENTROID  DATA  HI 
’/.================================= 

“/MATRIX  AND  VECTOR  CENTROID  METHOD 
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[Centroidl, Centroid2,Tl,T2,PDAlc,PDA2c,dXl,dX2,Centroidlr, . . . 
Centroid2r,dXlc,dX2c,Centroidlp,Centroid2p,Tlp,T2p]  =  . . . 

centroid(Max_Signall,Max_Signal2,Xl,X2,PDAl,PDA2, . . . 
ml ,m2 ,nl ,n2 , Tms_l ,Tms_2 ,PDAluse ,PDA2use) ; 


%================================================================ 

HI  ERROR  ANALYSIS  (delta  centroid2)  III 
%================================================================ 

[PDAlt,PDA2t,Centroidlt,Centroid2t,Tlt,T2t,dClt,dC2t]  =  ... 

centroiderror (ml ,m2,nl ,n2 ,PDAlc ,PDA2c ,X1 ,X2 ,dXl , dX2 ,T1 ,T2, . . . 
Centroidl ,Centroid2, shotnum,Cx, Cy, magtime ,PDAluse,PDA2use, . . . 
dXlc , dX2c) ; 


’/.============================================================= 

HI  ABEL  INVERSION  HI 

%============================================================= 

7„pdai 

[Lshelll ,Rshelll]  =  abell (XI ,nonzerol ,Tlt ,Centroidlt ,PDAlt , . . . 
magtime, model, Cx, Cy,  ..  .  °/0Mode  data  for  subplot 

PDA1 , shotnum,Tms_l ,  .  .  .  /(Contour  data  for  subplot 

PDAluse,ml);  /(PDA  in  use 


%PDA2 

[Lshell2 , Rshell2]  =  abel2 (X2 , nonzero2 , T2t , Centroid2t , PDA2t , . . . 
magtime, model, Cx, Cy,  ..  .  °/0Mode  data  for  subplot 

PDA2,shotnum,Tms_2,  .  .  .  /(Contour  data  for  subplot 

PDA2use,m2);  °/0PDA  in  use 


%===================================================================== 

l  DATA  PLOTS  Till 

dataplots (XI ,X2, Centroidl ,Centroid2 ,Tms_l ,Tms_2 ,T1 ,T2,PDA1 ,PDA2 , . . . 

Max_Signal_valuel ,Max_Signal_value2,nonzerol ,nonzero2 , shotnum, . . . 
PDAluse,PDA2use,Cx,Cy, magtime, Tip, T2p, Tit, T2t,Centroidlp, . . . 
Centroid2p,Centroidlt,Centroid2t,dXlc,dX2c, . . . 
Centroidlr,Centroid2r) ; 


%========================== 

°/„L00PING  SHOT  NUMBER  ENDING 
l  end 
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A. 2  tankaxis.m 


’/.=================================================================== 

°/o 

70Adam  Mads  on 
°/0l/18/05 

°/o 

“/.ARRAY  AXIS  VALUES  CALL  FUNCTION:  tankaxis.m 

% 

“/.Description:  This  call  function  properly  sets  the  chord  locations 
“/„  of  the  two  photodiode  arrays  and  references  them  to  the  proper 
“/.  x-  and  y-axis  contours. 

“/. 

’/.================================================================== 

“/. 

“/.begin  program 

function  [axisl , axis2]  =  tankaxis (nonzerol ,nonzero2) 
format  short 

if  size (nonzerol)  ==  [0  0] 

nzl  =  1;  al  =  0;  “/."a"  is  a  multiplier  to  return  zero  in  X  axis 

else 

nzl  =  nonzerol;  al  =  1; 

end 

if  size (nonzero2)  ==  [0  0] 
nz2  =1;  a2  =  0; 

else 

nz2  =  nonzero2;  a2  =  1; 

end 


“/.“/„“/„  Y-axis  array  locations  (cm)  °/.0/.0/. 

“/.  NOTE:  USE  ON  SIDE  MOUNTED  PDA’S  WHERE  CHANNEL  1  IS  ON  -Y 
“/.  (floor)  SIDE 

yraw  =  [-2.342,  -2.184,  -2.027,  -1.869,  -1.712,  -1.554,  -1.397,... 

-1.240,  -1.082,  -0.925,  -0.767,  -0.610,  -0.452,  -0.295,... 
-0.137,  0.020,  0.178,  0.335,  0.493,  0.650,  0.808,... 

0.965,  1.123,  1.280,  1.438,  1.595,  1.753,  1.910,... 

2.068,  2.225,  2.383,  2.540]; 
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yl  =  al*yraw; 
y2  =  a2*yraw; 


'/.REORDER  OF  AXIS  VALUES  FOR  180  DEGREE  PDA  ROTATION  ON  TANK 
7.  NOTE:  USE  ON  SIDE  MOUNTED  PDA’S  WHERE  CHANNEL  1  IS  ON  +Y 
7.  (CEILING)  SIDE 

ylflip  =  fliplr(yl) ; 
y2flip  =  fliplr(y2); 


7.7.7.  X-axis  array  locations  (cm)  7.7.7. 

'/.NOTE:  CURRENTLY  USING  SAME  VALUES  FOR  X  AND  Y  AXIS  VIEWS 
7,  NOTE:  -X  DIRECTION  IS  TOWARDS  PUMPING  SYSTEM  ...  therefore 
7.  *USE  fliplr(y)  VALUES  IF  CH.  1  IS  ON  +X  SIDE 

7.  *USE  yraw  VALUES  IF  CH.l  IS  ON  -X  SIDE 

'/.USER  CHANGED  PARAMETERS  BASED  ON  THAT  DAY’S  PDA  ORIENTATION  AND 
7.  REMOVES  UNUSED  CHANNELS 

axisl  =  yl(nzl);  '/.Currently  PDA1  @  PO  side  ...  ch.l  (-Y) 

axis2  =  y2flip(nz2);  '/.Currently  PDA2  @  PO  top  ...  ch.l  (+X) 
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A. 3  magneticmode. m 


7.================================================================== 

°/o 

°/0Adam  Mads  on 
°/„l/09/05 

7. 

°/„MAGNETIC  MODE  CALL  FUNCTION:  magneticmode . m 

% 

"/[Description:  This  call  function  takes  azimuthal  magnetic  probe 
7.  data  from  experiment  and  determines  the  normalized  mode  data 

°/„  and  Fourier  component  center  of  current  data  for  the  x-  and  y- 

°/„  axis.  Various  plotting  routines  relating  to  different 
°/„  components  of  this  data  are  listed  at  the  end  of  the  code. 

I 

7.================================================================== 

l 

°/0begin  program 

function  [magtime, mode, model ,Cx,Cy, cutoff]  =  ... 
magneticmode (shotnum,probelocat ion) 

format  short 

7.================================================================== 

“/.USER  INPUTS  IF  RUN  AS  A  PROGRAM  AND  NOT  AS  A  CALL  FUNCTION 
7.  shotnum  =  41103039; 

°/0Probe  Location 
7.  probelocation  =  0; 

°/0Access  data  from  MDS 

7.  [a ,  status]  =  mdsopen  ( ’  landau .  aa .  Washington .  edu :  :  zapmain  ’  ,  shotnum)  ; 
[a, status]  =  mdsopen ( ’alfven . aa. Washington . edu: : zapmain’ , shotnum) ; 

°/016384  data  points  from  16kB  digitzer  setting 
time  =  zeros (1 , 16384) ; 
ml  =  zeros (1 , 16384) ; 
mO  =  zeros (1 , 16384) ; 
s  =  zeros (1 , 16384) ; 
c  =  zeros (1 , 16384) ; 


7oWall  radius  data 
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r_wall=mdsvalue ( 5 \r_wall ’ ) ; 

°/„Magnetic  data  for  positive  z-axis  placement 
ml  =  mdsvalue(  [,\m_l_pJ ,num2str (probelocation)] ) ’ ; 
mO  =  mdsvalue(  [’ViLCLp’  ,num2str (probelocation)] ) ’  ; 
magtime  =  mdsvalue( [’dim_of (\m_l_p’ ,num2str (probelocation) 
magtime  =  le6*magtime; 

°/oEmpirical  value  cutoff  at  .2 

cutoff  =  ones (1 , length (magtime) ) ;  cutoff  =  . 2*cutoff; 

cutoff  1  =  100* (r_wall/2) *cutof f ; 

for  i  =  1: length (mO) 
if  m0(i)  >  .01 

mO (i)  =  mO (i)  ; 
else  m0(i)  =  .01; 
end 

end 

/(Normalized  mode  data 

%  Can  plot  this  for  axisymmetric  plot  of  center  of  current 
mode  =  ml . /mO ; 

model  =  100* (r_wall/2) *mode ; 

%Sine  and  Cosine  Data  from  Tree 

s  =  mdsvalue(  [,\s_ml_pJ ,num2str (probelocation)] ) ’ ; 
c  =  mdsvalue(  [,\c_ml_pJ ,num2str (probelocation)] ) ’ ; 

%Sine  and  Cosine  data  corrected  for  angle  to  axial 

“/„  probe  array  location 

sine  =  sin(atan2(s,c)-157.5*pi/180) ; 

cosine  =  cos(atan2(s,c)-157.5*pi/180) ; 

“/(Displacement  along  x  and  y  planes  (cm) 

Cx  =  100* (r_wall/2) *mode . *cosine; 

Cy  =  100* (r_wall/2) *mode . *sine ; 

°/„Tank  circle 
th  =  (0: 127) /128*2*pi ; 
x  =  5*cos (th) ; 
y  =  5*sin(th) ; 
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A. 4  centroid,  m 


’/.=================================================================== 

°/o 

70Adam  Mads  on 
%l/05/05 

°/o 

“/.CENTROID  CALCULATION  CALL  FUNCTION:  centroid. m 

% 

“/.Description:  The  centroid  function  determines  the  centroid  of 
“/.  emissivity  of  the  PDA  chord  integrated  signal  using  either 
7„  the  vector  and  matrix  method  (currently  used  for  error  analysis) 
7„  or  the  trapezoid  rule.  Centroid  information  is  found  for  both 
7.  PDA  arrays  (designated  by  the  1  and  2  after  stored  variables)  . 

“/„  This  program  also  truncates  non-peaked  data  points  and  large 
“/„  offsets  of  the  centroid  to  any  one  side. 

“/. 

"/.=================================================================== 

“/. 

“/„  begin  function 

function  [Centroidl ,Centroid2 ,T1 ,T2 ,PDAlc ,PDA2c ,dXl ,dX2, . . . 

Centroidlr,Centroid2r,dXlc,dX2c,Centroidlp,Centroid2p, . . . 
Tlp,T2p] . . . 

=  centroid(Max_Signall ,Max_Signal2 ,X1 ,X2 ,PDA1 ,PDA2 , . . . 
ml ,m2 ,nl ,n2,Tms_l ,Tms_2 ,PDAluse ,PDA2use) ; 

°/0Call  function  if  Trapezoid  part  is  uncommented 
"/(function  [Centroidl, Tl,PDAlc,dXl,Centroid_Trap]  =  ... 

7.  centroid ( shotnum, Max_Signal  1 , XI , PDA1  ,ml  ,nl  ,Tms_l) 

format  short 


"/.=========================================================== 

“/.“/„“/„  CENTROID  CALCULATIONS  USING  VECTOR  AND  MATRIX  METHOD  “/„“/„“/„ 
"/„=========================================================== 

7.  Used  for  Error  Analysis 

°/0Emissivity  with  forced  zeros  on  added  outside  channels 
Etotl  =  zeros (ml+2 ,nl) ; 

Etotl ( 1 ,  :  )  =  0 ; 

Etotl (end, : )  =  0 ; 
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Etotl (2 : end-1 , : )  =  PDA1; 

Etot2  =  zeros (m2+2,n2) ; 

Etot2(l ,  :  )  =  0; 

Etot2 (end, : )  =  0 ; 

Etot2 (2 : end-1 , : )  =  PDA2; 

°/„delta  XI  array 

dXlcent  =  zeros (1  ,ml-l) ;  70delta  XI  array  of  inside  channels 

dXl  =  zeros  (1  ,ml+l)  ;  °/0total  delta  XI  array 

for  i  =  l:ml-l 

dXlcent (i)  =  XI (i+1) -XI (i) ; 

end 

if  PDAluse  ==  0 
dXl  =  1; 

else 

dXl(l)  =  abs (mean (dXlcent) ) ; 
dXl(end)  =  abs (mean(dXlcent) ) ; 
dXl(2:end-l)  =  abs (dXlcent) ; 

end 

°/„delta  X2  array 

dX2cent  =  zeros  (1  ,m2-l) ;  °/0delta  X2  array  of  inside  channels 

dX2  =  zeros  (1  ,m2+l)  ;  °/0total  delta  X2  array 

for  i  =  l:m2-l 

dX2cent (i)  =  X2 (i+1) -X2 (i) ; 

end 

if  PDA2use  ==  0 
dX2  =  1; 

else 

dX2(l)  =  abs (mean (dX2cent) ) ; 
dX2(end)  =  abs (mean(dX2cent) ) ; 
dX2(2:end-l)  =  abs (dX2cent) ; 

end 

0/„Centered  XI  position  calculations 
Xtotl  =  zeros (1 ,ml+2) ; 
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Xlc  =  zeros (1 ,ml+l) ; 

Xtotl (2 : end-1)  =  XI; 

°/0Centered  X2  position  calculations 
Xtot2  =  zeros (1 ,m2+2) ; 

X2c  =  zeros (1 ,m2+l) ; 

Xtot2 (2 : end-1)  =  X2; 


°/0NOTE:  The  first  and  last  values  of  Xtot  are  dependent  on  the 
°/0  orientation  of  the  PDA  on  the  machine.  Make  sure  that  they  are 
°/0  to  the  left  and  right  of  the  first  and  last  values  of  X 
7.  respectively.  This  assumes  that  the  PDA  straddles  the  zero  axis. 

°/„PDAl 

if  XI (1)  <=  0 

Xtotl  (1)  =  Xl(l)-dXl(l); 

else 

Xtotl  (1)  =  XI (l)+dXl (1) ; 

end 

if  XI (end)  >=  0 

Xtotl (end)  =  XI (end)+dXl (end) ; 

else 

Xtotl (end)  =  XI (end) -dXl (end) ; 

end 

0/oPDA2 

if  X2(l)  <=  0 

Xtot2 (1)  =  X2(l)-dX2(l); 

else 

Xtot2 (1)  =  X2(l)+dX2(l) ; 

end 

if  X2(end)  >=  0 

Xtot2(end)  =  X2(end)+dX2(end) ; 

else 

Xtot2(end)  =  X2 (end) -dX2 (end) ; 

end 


/(Centered  XI  array 
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for  i  =  l:ml+l 

Xlc(i)  =  (  Xtotl(i)+Xtotl(i+l)  )/2; 

end 

°/„Centered  X2  array 
for  i  =  1 : m2+l 

X2c(i)  =  (  Xtot2(i)+Xtot2(i+l)  )/2; 

end 

%Centered  Emissivity  signals 
Elc  =  zeros (ml+1 ,nl) ; 

E2c  =  zeros (m2+l ,n2) ; 

for  i  =  l:ml+l 
for  j  =  1 : nl 

Elc(i,j)  =  (  Etotl (i , j )+Etotl (i+1 , j )  )/2; 

end 

end 

if  PDAluse  ==  0 
Elc  =  1; 

end 

for  i  =  1 : m2+l 
for  j  =  1 :n2 

E2c(i,j)  =  (  Etot2 (i , j )+Etot2(i+l , j )  )/2; 

end 

end 

if  PDA2use  ==  0 
E2c  =  1; 

end 

°/„CENTR0ID  CALCULATION 
Centroidlr  =  zeros(l,nl); 

Centroidlr  =  (dXl . *Xlc*Elc) . /(dXl*Elc) ; 

Centroid2r  =  zeros(l,n2); 

Centroid2r  =  (dX2 . *X2c*E2c) . /(dX2*E2c) ; 

°/„average  signal  of  outside  channels 
if  PDAluse  ==  0 
Endavgl  =  0; 
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else 

Endavgl  =  zeros (1 ,nl) ; 

Endavgl  =  (abs (PDA1 (1 , : ) )+abs(PDAl (ml , : ) ) )/2 ; 

end 

if  PDA2use  ==  0 
Endavg2  =  0; 

else 

Endavg2  =  zeros (l,n2); 

Endavg2  =  (abs (PDA2 (1 , : ) )+abs(PDA2(m2 , : ) ) )/2 ; 

end 

/(Truncate  PDA1  input  values  to  remove  non-peaked  data  sets 
p  =  0;  /(dummy  variable 

if  PDAluse  ==  0 

j=i; 

else 

j=nl; 

end 

for  i  =  1 :  j 

if  Max_Signall (i)  >=  3*Endavgl(i) 

P  =  [p,i] ; 

end 

end 

Centroidlp=Centroidlr  (p  (2 :  end) ) ;  /(returns  peaked  centroid  data 

Tip  =  Tms_l  (p(2  :end) )  ;  °/0time  of  peaked  data 

PDAlp  =  PDA1  ( :  ,p(2 :  end) )  ;  "/(signals  at  peaked  data  points 

/(Truncate  PDA2  input  values  to  remove  non-peaked  data  sets 
p  =  0;  /(dummy  variable 

if  PDA2use  ==  0 

j=i; 

else 

j=n2; 

end 

for  i  =  1 :  j 

if  Max_Signal2 (i)  >=  3*Endavg2 (i) 

P  =  [p,i] ; 
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end 


end 

Centroid2p=Centroid2r (p (2 : end) ) ;  %returns  peaked  centroid  data 

T2p  =  Tms_2  (p(2 :  end) ) ;  °/0time  of  peaked  data 

PDA2p  =  PDA2( : ,p(2:end)) ;  %signals  at  peaked  data  points 

70Truncate  PDA1  input  values  to  remove  points  where  centroid  is 
%  all  the  way  to  the  left  or  right  of  the  chords 
p  =  0;  °/„  dummy  variable 

for  i  =  1 : length(Tlp) 

°/0NOTE:  Watch  inequality  change  if  XI (1)  and  XI (end) 

°/„  centered  around  zero 

if  XI C 1)  <=  0 

if  Centroidlp(i)  >  X 1 ( 1 )  &  Centroidlp(i)  <  Xl(end) 

P  =  [p,i] ; 

end 

elseif  XI C 1 )  >=  0 

if  Centroidlp(i)  <  X 1 ( 1 )  &  Centroidlp(i)  >  Xl(end) 

P  =  [p,i] ; 

end 

end 

end 

if  PDAluse  ==  0 

Centroidl  =0;  T1  =  0;  PDAlc  =0;  dXl  =  0; 

else 

Centroidl=Centroidlp(p(2 : end) ) ;  %returns  correct  centroid  data 
T1  =  Tlp(p(2 :  end) )  ;  °/0time  of  truncated  data 

PDAlc  =  PDAlp( :  ,p(2 :  end) ) ;  /(signals  at  correct  data  points 

end 

/(Truncate  PDA2  input  values  to  remove  points  where  centroid  is 
%  all  the  way  to  the  left  or  right  of  the  chords 
p  =  0;  /(dummy  variable 

for  i  =  1 : length(T2p) 

°/0NOTE:  Watch  inequality  change  if  X2(l)  and  X2(end) 

7.  centered  around  zero 

if  X2(l)  <=  0 

if  Centroid2p(i)  >  X2(l)  &  Centroid2p(i)  <  X2(end) 

P  =  [p,i] ; 


end 
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elseif  X2(l)  >=  0 

if  Centroid2p(i)  <  X2(l)  &  Centroid2p(i)  >  X2(end) 

P  =  [p,i]  ; 

end 

end 

end 

if  PDA2use  ==  0 

Centroid2  =0;  T2  =  0;  PDA2c  =  0;  dX2  =  0; 

else 

Centroid2=Centroid2p(p(2 : end) ) ;  %returns  correct  centroid  data 
T2  =  T2p(p(2 :  end) ) ;  °/0time  of  truncated  data 

PDA2c  =  PDA2p( : ,p(2 : end) ) ;  %signals  at  correct  data  points 

end 

°/00btain  correctly  sized  center  to  center  dX  for  centroiderror  .m  code 
dXlc  =  []  ; 
dX2c  =  []  ; 

if  PDAluse  ==  0 
dXlc  =  0; 

else 

for  i  =  1 : length(Xlc)-l 

dXlc(i)  =  Xlc(i+1)-Xlc(i); 

end 

dXlc  =  abs(dXlc) ; 

end 

if  PDA2use  ==  0 
dX2c  =  0; 

else 

for  i  =  1 : length(X2c)-l 

dX2c (i)  =  X2c(i+1)-X2c(i); 

end 

dX2c  =  abs(dX2c); 

end 


’/.======================================================= 

HI  CENTROID  CALCULATIONS  USING  TRAPEZOID  RULE  HI 
’/.======================================================= 

°/„NEED  TO  CHANGE  VARIABLE  NAMES  TO  MATCH  NEW  DESIGNATIONS 
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%  Centroid_Trap  =  zeros (l,nl); 

%  dX_Trap=  zeros (1 ,ml-l) ; 

%  Xl_2  =  zeros (1 ,ml-l) ; 

%  Z  =  zeros (1 ,ml) ; 

%  for  j  =  l:nl 
7,  Z  =  PDA1  ( :  ,  j )  ; 

7»  for  i  =  l:ml-l 

7,  dX_Trap(i)  =  XI (i+1) -XI (i)  ; 

7,  Xl_2  (i)  =  (XI  (i+l)+Xl  (i)  )/2; 

7,  end 

l 

7»  Trap_den  =  zeros  (1  ,ml-l)  ; 

7»  Trap_num  =  zeros  (1  ,ml-l)  ; 

7»  for  i  =  2: ml 

7»  Trap_den(i-1)  =  . 5*dX_Trap(i-l) *(Z(i-l)+Z(i) )  ; 

7,  Trap_num(i-1)  =  .5*dX_Trap(i-l)*(Z(i-l)+Z(i))*Xl_2(i-l)  ; 

7,  end 

7,  Trap_D  =  sum(Trap_den)  ; 

7,  Trap_N  =  sum(Trap_num)  ; 

7, 

7»  Centroid_Trap(l,  j)  =  Trap_N/Trap_D; 

7,  end 
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A. 5  centroiderror.m 


7.=================================================================== 

7. 

70Adam  Mads  on 
°/0l/13/05 

7. 

“/.CENTROID  STANDARD  DEVIATION  CALL  FUNCTION:  centroiderror.m 

l 

"/(Description:  The  centroid  error  call  function  determines  the 
7.  standard  deviation  of  the  truncated  centroid  points  determined 
°/„  in  the  centroid. m  call  function  and  passed  to  PDAcode.m.  It 
°/„  plots  errorbars  of  these  deviations  at  every  50th  location  along 
7.  the  centroid  vector  and  plots  them  versus  time. 

7. 

7.=================================================================== 

7. 

“/.begin  program 

function  [PDAlt ,PDA2t ,Centroidlt ,Centroid2t ,Tlt ,T2t ,dClt ,dC2t]  =  ... 

centroiderror (ml ,m2,nl ,n2 ,PDAlc ,PDA2c ,X1 ,X2 ,dXl ,dX2 ,T1 ,T2, . . . 
Centroidl ,Centroid2, shotnum, Cx, Cy, magtime , . . . 

PDAluse , PDA2use , dXlc , dX2c) 

format  short 

7.7.7.  pda  l  USE  7.7.7. 

if  PDAluse  ==  0 

Centroidlt  =  0;  Tit  =  0;  PDAlt  =  0;  dClt  =  0; 

else 

“/.Delta  Signal  vector 

dVlc  =  zeros (ml-2 , length(PDAlc) ) ; 

for  i  =  l:ml-2 

for  j  =  1 : length(PDAlc) 

dVlc(i,j)  =  PDAlc(i+l, j)-PDAlc(i, j) ; 

end 

end 

dVl  =  zeros (ml , length(PDAlc) ) ; 
dV 1 ( 1 , : )  =  PDAlc (1 ,  :) ; 
dVl (end, :)  =  PDAlc (ml, :) ; 
dVl (2: end-1, :)  =  dVlc; 
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Clden  =  dXlc*PDAlc; 

dClsquared  =  zeros (1 , length (PDAlc) ) ; 

for  i  =  1 : length(PDAlc) 

dClsquared(i)=  (1/ (Clden (i) . ~2) ) . *.  .  . 

(dXlc . *X1 -Centro idl (i) . *dXlc) . ~2*dVl ( : , i) . ~2; 

end 

dCl  =  sqrt (dClsquared) ; 

/(Truncate  to  between  20  and  100  mu*s  for  plotting  purposes 
p  =  0;  7„dummy  variable 

for  i  =  l:length(Tl) 

if  T1 (i) >=20  &  T1 (i) <=100 
P  =  [p,i] ; 

end 

end 

Centroidlt=Centroidl (p(2 : end) ) ; 

Tit  =  T1 (p(2 : end) ) ;  /(Time  values  for  truncated  vector 

PDAlt  =  PDAlc (: ,p(2 : end) ) ;  /(Signals  for  truncated  matrix 
dClt  =  dCl  (p(2 :  end) )  ;  °/0Delta  centroid  at  truncated  vector 

/(Truncation  of  errorbars  to  sample  every  50th  location 
°/„  from  20-100  mu*s 

Tlsample  =  zeros(l,floor(length(Tlt)/50)) ; 
dClsample  =  zeros(l,floor(length(dClt)/50)) ; 

Centroidlsample  =  zeros(l,floor(length(Centroidlt)/50)) ; 
for  i  =  1 : 1 : length(Tlsample) 

Centroidlsample(i)  =  Centroidlt (i*50) ; 

Tlsample(i)  =  Tlt(i*50); 
dClsample(i)  =  dClt(i*50); 

end 

end 

%7o0/o  PDA  2  USE  7o7o7o 
if  PDA2use  ==  0 

Centroid2t  =  0;  T2t  =  0;  PDA2t  =  0;  dC2t  =  0; 

else 

/(Delta  Signal  vector 

dV2c  =  zeros (m2-2, length (PDA2c) ) ; 

for  i  =  l:m2-2 

for  j  =  1 : length (PDA2c) 
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dV2c(i,j)  =  PDA2c(i+l, j)-PDA2c(i, j) ; 

end 

end 

dV2  =  zeros (m2 , length(PDA2c) ) ; 
dV2 (1 ,  :)  =  PDA2c ( 1 ,  :)  ; 
dV2 (end, :)  =  PDA2c(m2, :) ; 
dV2 (2: end-1, :)  =  dV2c; 

C2den  =  dX2c*PDA2c; 

dC2squared  =  zeros (1, length (PDA2c) ) ; 

for  i  =  1 : length (PDA2c) 

dC2squared(i)=  (1/ (C2den(i) . ~2) ) . * . . . 

(dX2c . *X2-Centroid2 (i) . *dX2c) . ~2*dV2( : ,i) ,~2; 

end 

dC2  =  sqrt (dC2squared) ; 

/(Truncate  to  between  20  and  100  mu*s  for  plotting  purposes 
p  =  0;  °/„dummy  variable 

for  i  =  l:length(T2) 

if  T2 (i) >=20  &  T2 (i) <=100 
P  =  [p,i]  ; 

end 

end 

Centroid2t=Centroid2 (p (2 : end) ) ; 

T2t  =  T2 (p(2 : end) ) ;  %Time  values  for  truncated  vector 

PDA2t  =  PDA2c  ( :  ,p(2  :  end) )  ;  /(Signals  for  truncated  matrix 
dC2t  =  dC2  (p(2 :  end) )  ;  "/(Delta  centroid  at  truncated  vector 

"/(Truncation  of  errorbars  to  sample  every  50th  location 
°/„  from  20-100  mu*s 

T2sample  =  zeros(l,floor(length(T2t)/50)) ; 
dC2sample  =  zeros(l,floor(length(dC2t)/50)) ; 

Centroid2sample  =  zeros(l,floor(length(Centroid2t)/50)) ; 
for  i  =  1 : 1 : length(T2sample) 

Centroid2sample (i)  =  Centroid2t (i*50) ; 

T2sample(i)  =  T2t(i*50); 
dC2sample(i)  =  dC2t(i*50); 

end 


end 
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A. 6  abell.m 

%================================================================== 

% 

7,  A  dam  Mads  on 

‘/.1/05/05 

1 

°/„ABEL  INVERSION  CALL  FUNCTION:  abell.m 

1 

°/„N0TE:  Requires  centroid. m  and  centroiderror .m  to  run.  Also  uses 
*/.  magneticmode  .m  data  for  some  plots. 

% 

7oDescription:  The  Abel  routine  determines  the  radial  plasma  profile 
7»  from  the  chord  integrated  data  using  the  Abel  inversion  method. 
7»  The  abell.m  and  abel2.m  routines  are  identical  except  that  they 
7»  use  PDA1  or  PDA2  data.  Movie  and  plotting  routines  are  listed 
7»  at  the  end  of  this  code. 

7, 

%================================================================== 

7, 

7»begin  function 

7»  function  [Lshelll  ,Rshelll]  =  abell (XI ,nonzerol , Tit ,Centroidlt ,.. . 
7»  PDAlt,  .  .  .  7»Chord  data  used 

7»  magtime, model, Cx, Cy,  ..  .  7»Mode  data  for  subplot 

7»  PDA1 , shotnum,Tms_l ,  .  .  .  7«Contour  data  for  subplot 

7»  PDAluse,ml);  7oPDA  in  use 

7»f unction  without  magnetic  data 

7»  function  [Lshelll  ,Rshelll]  =  abell  (XI  ,nonzerol , Tit  ,Centroidlt ,..  . 
7»  PDAlt,  .  .  .  7»Chord  data  used 

7»  PDA1 , shotnum,Tms_l ,  .  .  .  7«Contour  data  for  subplot 

7»  PDAluse)  ;  7»PDA  in  use 

format  short 

7.BREAK  IF  PDA  NOT  IN  USE 
for  alpha  =  1 

if  PDAluse  ==  0 

Lshelll=0;  Rshelll=0; 
break 


else 
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°/0vector  and  array  creation 
n  =  length(Tlt) ; 
leftindex  =  []  ; 
rightindex  =  []  ; 
left  =  []  ; 
right  =  []  ; 

Lvect  =  []  ; 

Rvect  =  []  ; 
dleft  =  []  ; 
dright  =  []  ; 

Lshelll  =  zeros (length(nonzerol) ,n) ; 

Rshelll  =  zeros (length(nonzerol) ,n) ; 

°/0Ascending  order  check 
if  XI (1)  <  XI (2) 

XI  =  XI; 

PDAlt  =  PDAlt ; 
elseif  XI (1)  >  Xl(2) 

XI  =  fliplr(Xl); 

PDAlt  =  f lipud(PDAlt) ; 

end 

Zcreate  length  vector  to  define  shells  by  distance  from  centroid 
for  i  =  l:n 

signal  =  PDAlt (:,i); 

°/0left  side 

leftindex  =  find(Xl  <=  Centroidlt (i) ) ; 
left  =  XI (leftindex) ; 

Lvect  =  abs (lef t-Centroidlt (i) ) ; %chord  distances  from  centroid 
for  p  =  1 : length(Lvect)-l 

dleft (p)  =  abs (Lvect (p+1) -Lvect (p) ) ; 

end 

dleftavg  =  mean(dleft) ; 

°/0shell  radii  from  centroid  and  1st  chord  distance  from  cent, 
leftvect  =  [Lvect (D+dleftavg,  Lvect]; 

leftsignal  =  signal (leftindex) ;  %signal  at  each  chord 

LI  =  zeros (length(leftvect)-l) ; 

L2  =  zeros (length(leftvect)-l) ; 
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for  j  =  1 : length(lef tvect) -1 

for  k  =  1 : length(leftvect)-l 

if  leftvect (k) ~2-leftvect ( j+1) ~2  >  0 

Ll(j,k)  =  sqrt (leftvect (k) ~2-leftvect (j+1) ~2) ; 

else 

Ll(j.k)  =  0; 

end 

end 

for  k  =  1 : length(lef tvect)-l 

if  leftvect (k) “2-leftvect (j ) ~2  >  0 

L2(j,k)  =  sqrt (leftvect (k+1) “2-leftvect (j+1) ~2) ; 

else 

L2(j,k)  =  0; 

end 

end 

L  =  2*(L1-L2);  %left  side  length  matrix 

end 

LS  =  L\left signal; 

°/0Left  side  shell  emissivities  at  given  time  steps.  Row  order 
7.  from  outside  to  nearest  to  centroid 
Lshelll (1 : length(LS) , i)  =  LS; 

70right  side 

rightindex  =  find(Xl  >=  Centroidlt (i) ) ; 
right  =  XI (rightindex) ; 

Rvect  =  abs (right-Centroidlt (i) ) ;%chord  distances  from  centroid 
for  p  =  1 : length(Rvect) -1 

dright(p)  =  abs (Rvect (p+l)-Rvect (p) ) ; 

end 

drightavg  =  mean(dright) ; 

7„shell  radii  from  centroid  and  1st  chord  distance  from  cent, 
rightvect  =  [Rvect,  Rvect (end)+drightavg] ; 

rightsignal  =  signal (rightindex) ;  %signal  at  each  chord 

R1  =  zeros (length(rightvect)-l) ; 

R2  =  zeros (length(rightvect)-l) ; 

for  j  =  1 : length(rightvect)-l 

for  k  =  1 : length(rightvect)-! 
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if  rightvect (k+1) ~2-rightvect ( j ) ~2  >  0 

Rl(j,k)  =  sqrt (rightvect (k+1) ~2-rightvect(j)~2) ; 

else 

Rl(j,k)  =  0; 

end 

end 

for  k  =  1 : length (rightvect )-l 

if  rightvect (k) ~2-rightvect (j ) ~2  >  0 

R2(j,k)  =  sqrt (rightvect (k) ~2-rightvect(j) ~2) ; 

else 

R2 ( j , k)  =  0; 

end 

end 

R  =  2*(R1-R2);  °/0right  side  length  matrix 

end 

RS  =  R\right signal; 

°/0Right  side  shell  emissivities  at  given  time  steps.  Row  order 

7.  from  outside  to  nearest  to  centroid 

Rshelll (length (nonzerol)-length(RS)+l : end, i)  =  RS; 


end 

‘/.PDA  USE  LOOP 
end 


end 
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Appendix  B 

MATLAB  2-D  FOURIER  DECONVOLUTION  CODES 

Appendix  B  includes  the  MATLAB  code  and  calling  functions  dealing  with  two 
dimensional  Fourier  deconvolution.  Explanations  of  each  code  are  listed  in  their  ap¬ 
propriate  subsection.  Codes  with  movie  or  plotting  routines  are  not  included  in  these 
listings.  Shown  below  is  the  function  calling  sequence.  The  overall  code  FULL2d.m 
takes  chord  integrated  data  and  time  variables  from  the  PDAcode.m  function  and 
determines  the  deconvolved  radial  profile. 


CALLING  SEQUENCE 

•  FULL2d.m 


—  TDerror.m 


bestfit.m 
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B.l  FULL2d.m 


7. 

7.=================================================================== 

7. 

70Adam  Mads  on 

7.2/07/05 

7. 

°/„2D  FOURIER  DECONVOLUTION  FUNCTION:  FULL2d.m 

7. 

"/(Description:  FULL2d  is  the  overall  algorithm  for  the  2-D  Fourier 
°/0  deconvolution  of  PDA  emission.  Chord  integrated  data  or  radial 
°/0  synthetic  profiles  for  both  PDAs  are  called  at  the  beginning  of 
°/0  this  function.  Indeces  of  the  two  PDAs  are  matched  before 
°/0  generating  the  path  lengths  matricies  for  each  time  step.  At 
°/0  this  point,  TDerror.m  is  called  to  minimize  error  by  finding 

7.  the  chi~2  value.  The  final  offset  at  this  data  point  is  fed 

°/„  to  the  bestfit.m  routine  which  returns  the  radial  profile  at 
7.  that  time.  The  program  then  calls  the  next  set  of  chord  data 

7.  for  the  next  time  step.  Once  this  looping  is  completed,  the 

7.  plasma  radial  data  is  stored  in  a  multidimensional  array  and  the 

7.  offset  locations  in  vectors.  The  FULL2d.m  program  then 

7.  transforms  this  profile  from  the  plasma  to  the  lab  plane  and 

7.  plots  the  data. 

7. 

°/0NOTE:  This  program  takes  several  minutes  to  run.  Re-routing  the 
7.  call  functions  so  that  TDerror.m  calls  bestfit.m  for  the 
7.  majority  of  its  calculations  instead  of  going  through  its  own 
7.  looping  routine  might  improve  calculation  speed. 

% 

7.=================================================================== 

7. 


7obegin  program 
7o  clear  all;  close  all; 
format  short; 

7.  SHOT  VALUES 

7o  FEED  IN  XI  ,X2  ,PDAlt  ,PDA2t  ,Tlt  ,T2t ,  Centro idlt ,  Centroid2t  ,PDA1 ,  PDA2 
7o  USE  FULL  PDA  VALUES 
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%  %  CLEARING  UNNEEDED  VARIABLES 

clear  Centroidl  Centroid2  Centroidlp  Centroid2p  Centroidlr  . . . 
Centroid2r  Max_Signall  Max_Signal2  Max_Signal_valuel  .  .  . 
Max_Signal_value2  PDAlc  PDA2c  Tip  T2p  T1  T2  cutoff  . . . 
dClt  dC2t  dXl  dXlc  dX2  dX2c  mode  model 

warning  ( ;  off  5  ,  ’all’ )  °/0suppress  rank  deficient  warnings  that 

%  slow  down  the  code 


7,  7.BIT  NOISE  ERROR 

Eml  =  sum (PDA 1 (:, 1:100 ),2)/100; 

Em2  =  sum(PDA2 (:, 1:100 ),2)/100; 


%==================================== 

7» INDEX  MATCHING  OF  PDA1  AND  PDA2  DATA 
m  =  length(Tlt) ; 
n  =  length(T2t); 
storel  =  []  ; 
store2  =  []  ; 


if  m  <=  n 

p  =  1 ;  7odummy  variable 

for  i  =  1 : length(Tlt) 

Index  =  find(Tlt(i)  ==  T2t(l:end)); 
if  length(Index)  ~=  0 
storel(l,p)  =  i; 
store2(l,p)  =  Index; 
p  =  P+i; 

end 

end 

Tltt  =  Tlt(storel); 

PDAltt  =  PDAlt (:, storel) ; 

Centroidltt  =  Centroidlt (storel) ; 

T2tt  =  T2t(store2); 

PDA2tt  =  PDA2t ( : , store2) ; 

Centroid2tt  =  Centroid2t (store2) ; 
elseif  m  >  n 

p  =  1 ;  7odummy  variable 

for  i  =  1 : length(T2t) 

Index  =  find(T2t(i)  ==  Tit (1: end)); 
if  length(Index)  ~=  0 
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storel(l,p)  =  Index; 
store2(l,p)  =  i; 
p  =  P+i; 

end 

end 

Tltt  =  Tlt(storel); 

PDAltt  =  PDAlt ( : , storel) ; 
Centroidltt  =  Centroidlt (storel) ; 
T2tt  =  T2t(store2); 

PDA2tt  =  PDA2t ( : , store2) ; 
Centroid2tt  =  Centroid2t (store2) ; 

end 


7.============================================================= 

7.7.7.  PATH  LENGTH  MATRIX  GENERATION  7.7.7. 
%============================================================= 

°/0X-Axis  (Ex)  viewed  by  PDA  on  top  or  bottom  view  ports 
°/0Y-Axis  (Ey)  viewed  by  PDA  on  side  view  ports 

7. 

°/„N0TE:  WANT  ORIENTATION  FOR  PDA  PATH  LENGTHS  TO  BE  ASCENDING 
%  ORDER  IN  BOTH  TOP  AND  SIDE  VIEWS 

l 

°/„N0TE:  WANT  ORIENTATION  FOR  GRID  1  TO  CORRESPOND  TO  CHORD  1  OF 
%  PDA1  AND  PDA2  TO  BE  LOCATED  IN  -X,+Y  QUADRANT.  THEREFORE, 
%  TOP  PDA  NEEDS  CHORD  LOCATIONS  TO  GO  IN  ASCENDING  ORDER 
7.  (-X  ->  +X)  WHILE  SIDE  PDA  NEEDS  CHORDS  TO  GO  IN  DESCENDING 

7.  ORDER  (+Y  ->  -Y) 

°/0PDA  LOCATION 

7.  PDA1 :  XI  :P0-side 

7.  PDA2:  X2  :PO-vertical 

if  XI  (1)  >  XI  (2)  °/0side  view  (ascending  check) 

XI  =  fliplr(Xl); 

PDAltt  =  f lipud(PDAltt) ; 

end 

if  X2(l)  >  X2(2)  7otop  view  (ascending  check) 

X2  =  f liplr (X2) ; 

PDA2tt  =  f lipud(PDA2tt) ; 


end 
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’/.USER  INPUT 

Ex  =  PDA2tt;  °/„Ex  =  signal  from  TOP  VIEW 

Ax  =  X2; 

Centx  =  Centroid2tt; 

Ey  =  PDAltt;  °/„Ey  =  signal  from  SIDE  VIEW 

Ay  =  XI; 

Centy  =  Centroidltt; 


’/.================================================================ 

°/„CENTER  OF  PROFILE  FITTING  ROUTINE 
’/.================================================================ 

er_array  =  zeros ( 15, length(Tltt) ) ; 

Xoffset  =  zeros (1 , length(Tltt) ) ; 

Yoffset  =  zeros (1 , length(Tltt) ) ; 
r_array  =  zeros (15 , length (Tltt) ) ; 

for  i  =  1 : length(Tltt) 

“/.clearing  looping  variables 
clear  offsetX  offsetY  ErrorMesh 

Xguess  =  Centx(i);  “/.Initial  Profile  Center  guess 

Yguess  =  Centy(i); 

deltaX  =  1.0; 

deltaY  =1.0; 

“/.USE  FOR  SYNTHETIC  DATA 
7.  deltaX  =  .5; 

7.  deltaY  =  .5; 

dX  =  zeros (1 , length(Ax) -1) ; 
dY  =  zeros (1 , length(Ay) -1) ; 

for  j  =  1 : length(Ax) -1 

dX(j)  =  Ax(j+1)-Ax(j)  ; 

end 

for  j  =  1 : length(Ay) -1 

dY(j)  =  Ay ( j+1) -Ay ( j ) ; 

end 

dXavg  =  mean(dX) ; 
dYavg  =  mean(dY) ; 
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Chordavg  =  (dXavg+dYavg) /2 ; 

deltaX  =  . 25*length(Ax) *dXavg; 
deltaY  =  . 25*length(Ay) *dYavg; 


°/„LENGTH  MATRIX  GENERATION 
signal 1  =  Ex( : , i)  ; 
signal2  =  Ey( : ,i)  ; 
signal  =  [signall;  signal2] ; 
offsetX  =  zeros(9,l); 
offsetY  =  zeros(9,l); 

MinError  =  []  ; 

ERR1  =  []  ; 

ERR2  =  []  ; 

delta  =  max([deltaX  deltaY]); 

j  =  i; 

while  delta  >  . 05 

"/_========================================================= 

°/„MINIMUM  ERROR  LOOP  IN  X  AND  Y  (DEFINED  BY  DELTA) 
/„========================================================== 

[dEm_s,Xguessnew,Yguessnew,deltaXnew,deltaYnew,min_error, . . . 
of f setl , of f set2 , Error ,errl , err2 ,Ecal2]  =  ... 
TDerror(Ax, Ay, Xguess,Yguess,deltaX, deltaY, . . . 
signal ,Eml ,Em2) ; 

MinError(j)  =  min_error; 

ERRl(j)  =  errl; 

ERR2(j)  =  err2; 
for  k  =  1:9 

offsetX(k,j)  =  offsetl(k); 
offsetY(k,j)  =  offset2(k); 

ErrorMesh(k, j )  =  Error (k); 

end 

Xguess  =  Xguessnew; 

Yguess  =  Yguessnew; 

delta(j)  =  max([deltaX  deltaY]); 


deltaX  =  deltaXnew; 
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deltaY  =  deltaYnew; 


if  j  >=  500 
break 

end 

j  =  j+i; 

end 

finalerror  =  min(MinError) ; 

P  =  min (find (min (min (ErrorMesh) )  ==  ErrorMesh)); 

Xoff  =  offsetX(P); 

Yoff  =  offsetY(P); 

?.============================================================== 

“/RETURN  ARRAY  OF  BEST  RADIAL  PROFILES 
’/.============================================================== 

[erT,M,Ml ,M2,M3,r,X,Y,erraw]  =  bestf it (Ax, Ay , Xoff , Yoff , signal) ; 

%TIME  STEP  SAVES 
Xoffset(i)  =  Xoff; 

Yoffset(i)  =  Yoff; 

er_array (1 : length(erT) , i)  =  erT; 

r_array (1 : length (r) , i)  =  r; 

end 


“/=============================================================== 

“/.CONVERSION  TO  POLAR  AND  CARTESIAN  ARRAYS 
“/=============================================================== 

er_arr  =  []  ; 

er_arr  =  [er_array(l , : ) ; zeros(l , length (er_array) ) ; . . . 
zeros (1 , length (er_array) ) ; er_array (2 : end, : )] ; 

[m  n]  =  size (er_arr) ; 

q  =  (0 :  pi/16 :  2*pi) ;  “/.Number  of  Theta  points 

maxZ  =  zeros (1 , length(Tltt) ) ; 

“/Cartesian  Mesh  with  Center  radii  array  generation 
x  =  zeros (m/3+1 , length(q) , length(Tltt) ) ; 
y  =  zeros (m/3+1 , length(q) , length(Tltt) ) ; 
z  =  zeros (m/3+1 , length(q) , length(Tltt) ) ; 
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xn  =  zeros (m/3+1 , length(q) , length(Tltt) ) ; 
yn  =  zeros (m/3+1 , length(q) , length(Tltt) ) ; 
zn  =  zeros (m/3+1 , length(q) , length(Tltt) ) ; 

for  k  =  1 : length(Tltt) 

°/0Polar  Coordinate  Values 
Rho  =  zeros (m/3,  length (q) ) ; 

Theta  =  zeros (m/3, length (q) ) ; 

Z  =  zeros (m/3 , length(q) ) ; 

p  =  1; 

for  i  =  1 : length (r_array (f ind(r_array( : ,k) ) ,k) ) 
for  j  =  l:length(q) 

Theta(i, j)  =  q(j) ; 

Rho(i,j)  =  r_array (i ,k) ; 

Z(i,j)  =  er_arr (p,k)+er_arr (p+1 ,k)*cos(q(j))+. . . 
er_arr(p+2,k)*sin(q(j)) ; 

end 

P  =  P+3; 

end 

delta_r_ray=. . . 

(r_array (length (find (r_array ( : ,k) ) ) ,k)-r_array (1 ,k) )/ . . . 

(r_array (length (find (r_array( : ,k) ) ) ,k) -1) ; 

for  i  =  length(r_array(find(r_array( : ,k)) ,k))+l :m/3 
for  j  =  l:length(q) 

Theta(i, j)  =  q(j) ; 

Rho(i,j)  =  r_array(i-l ,k)+delta_r_ray; 

end 

P  =  P+3; 

end 

meanZ(k)  =  mean(Z(l ,  : ) ) ; 

“/CARTESIAN  COORDINATE  TRANSFORMATION 
°/„  only  needed  if  Theta  and  Rho  are  vectors 

°/„  [th_bogus ,r_bogus]  =  meshgrid(Theta,Rho) ; 

°/„  [xc,yc]  =  pol2cart(th_bogus,r_bogus)  ; 


[xc( : , : ,k) ,yc( : , : ,k) ,zc( : , : ,k)]  =  pol2cart (Theta, Rho ,Z) ; 
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°/0Adding  position  at  center 

xl  =  zeros (1 , length(q) )  ; 

yl  =  zeros (1 , length(q) ) ; 

zl  =  ones (1 , length (q) ) *meanZ(k) ; 

x(l, : ,k)  =  xl ; 

y(l. : >k)  =  yl; 

z(l , : ,k)  =  zl ; 
x(2:end,:,k)  =  xc(:,:,k); 
y (2 : end, : ,k)  =  yc(: , : ,k) ; 
z(2:end,:,k)  =  zc(:,:,k); 

xn(:,:,k)  =  x( : , : ,k)+Xoffset(k) ; 
yn(:,:,k)  =  y( : , : ,k)+Yoff set(k) ; 
zn( : , : ,k)  =  z( : , : ,k) ; 

end 
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B.2  TDerror.m 


’/.=================================================================== 

°/o 

°/0Adam  Mads  on 

%2/15/05 

*/. 

°/„TW0  DIMENSIONAL  FOURIER  DECONVOLUTION  ERROR  FINDER:  TDerror.m 

% 

"/(Description:  TDerror.m  is  a  calling  function  of  TDerror.m.  Given 
7.  an  initial  guessing  range,  it  will  find  the  minimum  error  in 
7„  terms  of  the  chi~2  values  of  a  given  profile  to  minimize  error 

7„  and  find  the  best  shell  profile  center.  Data  points  are  picked 

°/„  from  a  9  unit  grid  using  the  current  guess  value  and  distances 
7.  between  guesses.  If  the  minimum  returned  value  of  the  error  is 
°/„  at  the  center,  the  routine  lowers  the  spacing  of  the  guess.  If 
°/„  it  is  to  the  left  or  right,  it  moves  the  Xguess  or  Yguess  value 
7.  in  that  direction. 

% 

/„=================================================================== 

function  [dEm_s,Xguessnew,Yguessnew,deltaXnew,deltaYnew,min_error, . . . 
of fsetl , off set2, Error , errl , err2 ,Ecal2]  =  ... 

TDerror (Ax, Ay , Xguess , Yguess ,deltaX, deltaY, signal ,Eml ,Em2) ; 

format  short 

Error  =  zeros(3,3); 
errl  =  []  ; 
err2  =  []  ; 
min_error  =  []  ; 

for  g  =  1:9 

offsetl  =  [Xguess-deltaX, 

Xguess-deltaX, 

Xguess-deltaX, 

offset2  =  [Yguess+deltaY, 

Yguess , 

Yguess-deltaY, 


Xguess , 

Xguess , 

Xguess , 

Yguess+deltaY, 
Yguess , 

Yguess-deltaY, 


Xguess+deltaX; . . . 
Xguess+deltaX;  .  .  . 
Xguess+deltaX]  ; 

Yguess+deltaY;  .  .  . 
Yguess ; . . . 
Yguess-deltaY] ; 


er  =  [] 
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erT  =  []  ; 

°/„X-AXIS  CHORDS  AND  Y-AXIS  CHORDS 

“/oNOTE:  THESE  CHORD  LOCATIONS  SHOULD  BE  USED  FOR  REGULATING 
7.  SHELL  SIZES 

X  =  Ax-off setl(g) ; 

Y  =  Ay-off set2(g) ; 

°/„DELTA-X  and  DELTA-Y  vectors 
dX  =  zeros (1 , length(X)-l) ; 
dY  =  zeros (1 , length(Y)-l) ; 

for  j  =  1 : length(X)-l 
dX( j )  =  X( j+1) -X(j) ; 

end 

for  j  =  1 : length(Y)-l 
dY ( j )  =  Y ( j+1) -Y (j) ; 

end 

dXavg  =  mean(dX) ; 
dYavg  =  mean(dY) ; 

Chordavg  =  (dXavg+dYavg) /2 ; 

°/„SHELL  CREATION  (r) 

/(Minimum  set  by  empirical  distance  of  .01  relative  from  chords 
Xmin  =  find(abs(X)  >  .01); 

Xmin  =  X(Xmin); 

Ymin  =  find(abs(Y)  >  .01); 

Ymin  =  Y (Ymin) ; 

minimum  =  [abs(Xmin)  abs(Ymin)] ; 
innershell  =  Chordavg  +  min (minimum) ; 

/(Maximum  set  by  2nd  smallest  X  Y  Chord  value  greater  than  .  1 
maximum  =  sort ( [abs (Xmin (1) )  abs (Xmin (end) )  ... 

abs(Ymin(l))  abs (Ymin (end) )]) ; 
outershell  =  Chordavg  +  maximum(2) ; 

r  =  [innershell :  Chordavg :  outershell]  ;  /(Shell  radii 

Radii  =  zeros (1 , length(r)+l) ; 

Radii  =  [0,r] ; 
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T1  =  zeros (length(X) , length(r) ) ; 

T2  =  zeros (length(X) , length(r) ) ; 

51  =  zeros (length (Y) , length(r) ) ; 

52  =  zeros (length (Y) , length(r) ) ; 

0/oTOP/BOTTOM  VIEW 
for  j  =  l:length(X) 

for  k  =  l:length(r) 

if  Radii (k+1) ~2-X(j ) ~2  >  0 

Tl(j,k)  =  sqrt (Radii (k+1) ~2-X(j ) ~2) ; 

else 

T1 (j ,k)  =  0; 

end 

end 

for  k  =  l:length(r) 

if  Radii (k) ~2-X(j ) ~2  >  0 

T2(j,k)  =  sqrt (Radii (k) ~2-X(j ) ~2) ; 

else 

T2 ( j , k)  =  0; 

end 

end 

T  =  2* (T1-T2) ; 

end 

°/„SIDE/SIDE  VIEW 
for  j  =  l:length(Y) 

for  k  =  l:length(r) 

if  Radii (k+1) ~2-Y(j)~2  >  0 

Sl(j,k)  =  sqrt (Radii (k+1) ~2-Y(j ) ~2) ; 

else 

Sl(j,k)  =  0; 

end 

end 

for  k  =  l:length(r) 

if  Radii (k) ~2-Y(j ) ~2  >  0 

S2(j,k)  =  sqrt (Radii (k) ~2-Y(j ) ~2) ; 

else 

S2(j,k)  =  0; 

end 

end 

S  =  2*(S1-S2); 
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end 

“/oMIDPOINT  GENERATION 

MIDt_pos  =  (Tl+T2)/2;  %Y  value  midpoints  for  +Y  quadrants  II, I 

MIDs_neg  =  (Sl+S2)/2;  °/„X  value  midpoints  for  -X  quadrants  III, II 

MIDt_neg  =  -MIDt_pos;  %Y  value  midpoints  for  -Y  quadrants  III, IV 

MIDs_pos  =  -MIDs_neg;  °/„X  value  midpoints  for  +X  quadrants  IV, I 

Xn  =  find(X  <  0)  ;  °/0index  of  negative  X  values 

Xp  =  find(X  >=  0)  ;  °/0index  of  positive  X  values 

Yn  =  find(Y  <  0)  ;  °/0index  of  negative  Y  values 

Yp  =  find(Y  >=  0)  ;  °/0index  of  positive  Y  values 

“/oTHETA  VALUES  BY  QUADRANT  FOR  TOP  AND  SIDE  VIEWS 
7.  These  values  are  represented  in  column  order  according  to 
7„  the  corresponding  values  of  the  X  and  Y  vectors  respectively 
7„  EXAMPLE:  Any  Theta  values  corresponding  to  X(l)  would  be  in 

°/„  quadrant  II  or  III  due  to  ascending  order  of  X  vector  and 

°/„  would  be  in  the  first  row  of  the  Theta  matrix  for  those 

°/„  quadrants.  The  midpoint  values  used  would  relate  to  the 

°/„  quadrant  sign  (+/-)  . 

I 

°/0Top  view  quadrant  I  is  upper  right  corner 

Theta_tl  =  [] ;  %Top  view  Theta  values  for  quadrant  I 

Theta_tII  =  [] ;  %Top  view  Theta  values  for  quadrant  II 

Theta_tIII  =  [] ;  %Top  view  Theta  values  for  quadrant  III 

Theta_tIV  =  [] ;  %Top  view  Theta  values  for  quadrant  IV 

°/0Side  view  quadrant  I  is  upper  left  corner  due  to  90  deg 
7.  theta  rotation 

Theta_sl  =  []  ;  °/0Side  view  Theta  values  for  quadrant  I 

Theta_sII  =  []  ;  °/0Side  view  Theta  values  for  quadrant  II 

Theta_sIII  =  []  ;  °/0Side  view  Theta  values  for  quadrant  III 

Theta_sIV  =  []  ;  °/0Side  view  Theta  values  for  quadrant  IV 

for  k  =  l:length(r) 

for  j  =  l:length(Xn) 

Theta_tII ( j ,k)  =  atan2(MIDt_pos(j ,k) ,X(j)) ; 

Theta_tIII (j ,k)  =  atan2(MIDt_neg(j ,k) ,X(j)) ; 

end 

for  j  =  l:length(Xp) 

Theta_tl ( j ,k)  =  atan2(MIDt_pos (length (Xn)+j ,k) , . . . 
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X ( length (Xn) + j ) ) ; 

Theta_tIV(j ,k)  =  atan2(MIDt_neg (length (Xn)+j ,k) , . . . 
X(length(Xn)+j) ) ; 

end 

for  j  =  l:length(Yn) 

Theta_sIV(j ,k)  =  atan2(Y(j) ,MIDs_pos(j ,k)) ; 
Theta_sIII ( j ,k)  =  atan2 (Y ( j ) ,MIDs_neg( j ,k) ) ; 

end 

for  j  =  l:length(Yp) 

Theta_sl (j ,k)  =  atan2 (Y(length(Yn)+j ) , . . . 

MIDs_pos (length(Yn)+j ,k) ) ; 

Theta_sII (j ,k)  =  atan2(Y(length(Yn)+j) , . . . 
MIDs_neg(length(Yn)+j ,k) ) ; 

end 

end 


"/(Hemisphere  Theta  Matricies 
Theta_tN  =  [Theta_tII;  Theta_tl] ; 
Theta_tS  =  [Theta_tIII;  Theta_tIV] ; 
Theta_sE  =  [Theta_sIII;  Theta_sII] ; 
Theta_sW  =  [Theta_sIV;  Theta_sl] ; 


%Top  view  upper  half 
%Top  view  lower  half 
%Side  view  left  half 
%Side  view  right  half 


°/„SUBMATRIX  GENERATION 
°/0Radially  dependent  submatrix 
Ml  =  [T;S]  ; 

[a,b]  =  size (Ml) ; 


/(Cosine  dependent  submatrix 

topC  =  [cos (Theta_tN)+cos (Theta_tS) ]  ; 

sideC  =  [cos(Theta_sE)+cos(Theta_sW)] ; 

fullC  =  [topC; sideC] ; 

f or  j  =  1 : a 

for  k  =  l:length(r) 

M2 ( j , k)  =  (Ml ( j ,k) /2) . *f ullC ( j , k) ; 

end 

end 


°/0Sine  dependent  submatrix 
topS  =  [sin(Theta_tN)+sin(Theta_tS) ]  ; 
sideS  =  [sin(Theta_sE)+sin(Theta_sW)] ; 
fullS  =  [topS ; sideS] ; 
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for  j  =  1 :  a 

for  k  =  l:length(r) 

M3(j,k)  =  (Ml(j,k)/2).*fullS(j,k); 

end 

end 

“/.FULL  GEOMETRIC  MATRIX 
M  =  zeros (a,3*length(r) ) ; 

p  =  1; 

for  j  =  l:length(r) 

M(: ,p)  =  Ml ( : , j ) ; 

M(: ,p+l)  =  M2 ( : ,j) ; 

M( : ,p+2)  =  M3 ( : , j ) ; 

P  =  P+3; 

end 

"/.AXISYMMETRIC  1ST  SHELL 

M ( : ,2:3)  =  [] ;  "/.force  axisymmetric  first  shell 

[number_E,number_e] =size (M) ; 

°/0axi symmetric ,  cosine,  and  sine  component  indicies 

10  =  2:3:number_e;  °/0axisymmetric  indicies 

11  =  3:3:number_e;  °/0cosine  indicies 

12  =  4:3:number_e;  °/0sine  indicies 

10=  [1,10];  7„axisymmetric  for  truncated  M  matrix 

“/.LEAST  SQUARES  FIT 
erraw  =  zeros (number_e, 1) ; 

Ecal  =  zeros (number_E, 1) ; 

erraw  =  M\signal; 

“/.Constraint  on  Cosine  and  Sine  components 
eraxi  =  erraw(IO) ; 

ervar  =  2*sqrt (erraw(Il) . ~2+erraw(I2) . ~2) ; 
for  j  =  l:length(I2) 

if  eraxi (j+1)  <  ervar(j) 
eraxi(j+l)  =  ervar(j); 

end 

er(I0(j+l))  =  eraxi(j+l); 


end 
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er (1)  =  eraxi (1) ; 
er(Il)  =  erraw(Il); 
er(I2)  =  erraw(I2); 
erT  =  er’  ; 

°/0Error  calculation 

Ecal  =  M*erraw;  °/0Calculated  Signal 

Ecal2  =  M*erT;  °/0Calculated  Signal  with  Constraint 

°/0CHI~2  Error  Analysis 
"/Simulated  Data 
dEm_s  =  [Eml ; Em2] . ~2 ; 

errl(g)  =  sum( (Ecal-signal) . ~2 ,/dEm_s) ; 
err2(g)  =  sum( (Ecal2-signal) . ~2 . /dEm_s) ; 

7.  /EQUAL  WEIGHTING 

Error (g)  =  (sum( (Ecal-signal) . ~2 . /dEm_s)+ . . . 

sum( (Ecal2-signal) . "2 . /dEm_s) ) /2; 

end 

/Error  Minimization 

index_min  =  find(abs (Error)  ==  min (min (abs (Error) ))) ; 

/Reduce  Delta  Variation 
switch  index_min(l) 

/Shift  Xguess  and  Yguess 
case  1 

Xguessnew  =  (Xguess-deltaX)/2; 

Yguessnew  =  (Yguess+deltaY)/2; 
deltaXnew  =  deltaX*.9; 
deltaYnew  =  deltaY*.9; 
case  2 

Xguessnew  =  (Xguess-deltaX)/2 ; 

Yguessnew  =  Yguess; 
deltaXnew  =  deltaX*.9; 
deltaYnew  =  deltaY*.9; 
case  3 

Xguessnew  =  (Xguess-deltaX)/2; 

Yguessnew  =  (Yguess-deltaY)/2; 
deltaXnew  =  deltaX*.9; 
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deltaYnew  =  deltaY*.9; 
case  4 

Xguessnew  =  Xguess; 

Yguessnew  =  (Yguess+deltaY)/2; 
deltaXnew  =  deltaX*.9; 
deltaYnew  =  deltaY*.9; 

°/„Centered  DeltaX  and  DeltaY  reduction 
case  5 

Xguessnew  =  offsetl(5); 

Yguessnew  =  offset2(5); 
deltaXnew  =  deltaX*.7; 
deltaYnew  =  deltaY*.7; 

“/oShift  Xguess  and  Yguess 
case  6 

Xguessnew  =  Xguess; 

Yguessnew  =  (Yguess-deltaY) /2; 
deltaXnew  =  deltaX*.9; 
deltaYnew  =  deltaY*.9; 
case  7 

Xguessnew  =  (Xguess+deltaX)/2; 
Yguessnew  =  (Yguess+deltaY)/2; 
deltaXnew  =  deltaX*.9; 
deltaYnew  =  deltaY*.9; 
case  8 

Xguessnew  =  (Xguess+deltaX)/2; 
Yguessnew  =  Yguess; 
deltaXnew  =  deltaX*.9; 
deltaYnew  =  deltaY*.9; 
case  9 

Xguessnew  =  (Xguess+deltaX)/2; 
Yguessnew  =  (Yguess-deltaY) / 2 ; 
deltaXnew  =  deltaX*.9; 
deltaYnew  =  deltaY*.9; 
otherwise 

disp(; ERROR’) 

end 

min_error  =  min (min (abs (Error) )) ; 
errl  =  min(errl) ; 
err2  =  min(errl) ; 
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B.3  bestfit.m 


’/.=================================================================== 

°/o 

°/0Adam  Mads  on 

%2/15/05 

*/. 

°/„TW0  DIMENSIONAL  FOURIER  DECONVOLUTION  ERROR  FINDER:  TDerror.m 

% 

"/(Description:  Bestfit.m  is  a  calling  function  of  FULL2d.m.  It  takes 
7.  the  Xoff  and  Yoff  location  chosen  by  the  TDerror.m  function  as 

7„  the  best  point  to  deconvolve  the  profile.  It  then  creates  a 

7„  path  length  matrix  and  performs  a  least  squares  fit  on  the  chord 
°/„  integrated  data  to  produce  a  final  radial  profile  in  the  plasma 
7.  reference  frame.  This  data  is  then  fed  back  into  the  FULL2d.m 
°/„  code. 

°/o 

•/.=================================================================== 

function  [erT,M,Ml,M2,M3,r,X,Y,erraw]  =  ... 
bestf it (Ax, Ay , Xoff , Yoff , signal) ; 

format  short 

Xcenter  =  Xoff; 

Ycenter  =  Yoff ; 

°/„X-AXIS  CHORDS  AND  Y-AXIS  CHORDS 

°/„N0TE:  THESE  CHORD  LOCATIONS  SHOULD  BE  USED  FOR  REGULATING 
7.  SHELL  SIZES 

X  =  Ax-Xcenter; 

Y  =  Ay-Ycenter; 

°/„DELTA-X  and  DELTA-Y  vectors 
dX  =  zeros (1 , length(X)-l) ; 
dY  =  zeros (1 , length(Y)-l) ; 

for  j  =  1 : length(X)-l 
dX( j )  =  X( j+1) -X(j) ; 

end 

for  j  =  1 : length(Y)-l 
dY(j)  =  Y ( j+1) -Y (j)  ; 
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end 

dXavg  =  mean(dX) ; 
dYavg  =  mean(dY) ; 

Chordavg  =  (dXavg+dYavg) /2 ; 

“/.SHELL  CREATION  (r) 

“/.Minimum  set  by  an  empirical  distance  of  .  1  relative  from  any  chord 
Xmin  =  find(abs(X)  >  .01); 

Xmin  =  X(Xmin); 

Ymin  =  find(abs(Y)  >  .01); 

Ymin  =  Y (Ymin) ; 

minimum  =  [abs(Xmin)  abs(Ymin)] ; 
innershell  =  Chordavg;  +  min (minimum) ; 

“/.Maximum  set  by  2nd  smallest  X  Y  Chord  value  greater  than  .  1 
maximum  =  sort ( [abs (Xmin (1) )  abs (Xmin (end) )  ... 

abs(Ymin(l))  abs (Ymin (end) )]) ; 
outershell  =  Chordavg  +  maximum(2); 

r  =  [innershell :  Chordavg :  outershell]  ;  °/0Shell  radii 

Radii  =  zeros (1 , length(r)+l) ; 

Radii  =  [0,r] ; 

T1  =  zeros (length(X) , length(r) ) ; 

T2  =  zeros (length(X) , length(r) ) ; 

51  =  zeros (length (Y) , length(r) ) ; 

52  =  zeros (length (Y) , length(r) ) ; 

%T0P/B0TT0M  VIEW 
for  j  =  l:length(X) 

for  k  =  l:length(r) 

if  Radii (k+1) ~2-X(j ) "2  >  0 

Tl(j,k)  =  sqrt (Radii (k+1) ~2-X(j ) ~2) ; 

else 

T1 ( j , k)  =  0; 

end 

end 

for  k  =  l:length(r) 

if  Radii (k) ~2-X(j ) ~2  >  0 

T2(j,k)  =  sqrt (Radii (k) ~2-X(j ) ~2) ; 
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else 

T2(j,k)  =  0; 

end 

end 

T  =  2*(T1-T2) ; 

end 

7SIDE/SIDE  VIEW 
for  j  =  l:length(Y) 

for  k  =  l:length(r) 

if  Radii (k+1) ~2-Y( j ) ~2  >  0 

Sl(j,k)  =  sqrt (Radii (k+1) ~2-Y(j ) ~2) ; 

else 

Sl(j,k)  =  0; 

end 

end 

for  k  =  l:length(r) 

if  Radii (k) ~2-Y(j ) ~2  >  0 

S2(j,k)  =  sqrt (Radii (k) ~2-Y(j ) ~2) ; 

else 

S2(j,k)  =  0; 

end 

end 

S  =  2* (S1-S2) ; 

end 

7.MIDP0INT  GENERATION 

MIDt_pos  =  (Tl+T2)/2;  °/0Y  value  midpoints  for  +Y  quadrants  II, I 

MIDs_neg  =  (Sl+S2)/2;  °/0X  value  midpoints  for  -X  quadrants  III, II 

MIDt_neg  =  -MIDt_pos;  °/0Y  value  midpoints  for  -Y  quadrants  III, IV 

MIDs_pos  =  -MIDs_neg;  °/0X  value  midpoints  for  +X  quadrants  IV, I 

Xn  =  find(X  <  0)  ;  70index  of  negative  X  values 

Xp  =  find(X  >=  0)  ;  70index  of  positive  X  values 

Yn  =  find(Y  <  0)  ;  Tiindex  of  negative  Y  values 

Yp  =  find(Y  >=  0)  ;  Tiindex  of  positive  Y  values 

7.THETA  VALUES  BY  QUADRANT  FOR  TOP  AND  SIDE  VIEWS 
T  These  values  are  represented  in  column  order  according  to  the 
T  corresponding  values  of  the  X  and  Y  vectors  respectively 
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%  EXAMPLE.  Any  Theta  values  corresponding  to  X(l)  would  be  in 

%  quadrant  II  or  III  due  to  ascending  order  of  X  vector  and  would 

%  be  in  the  first  row  of  the  Theta  matrix  for  those  quadrants. 

1  The  midpoint  values  used  would  relate  to  the  quadrant  sign  (+/-) 

% 

%Top  view  quadrant  I  is  upper  right  corner 

Theta_tl  =  []  ;  7>Top  view  Theta  values  for  quadrant  I 

Theta_tII  =  []  ;  7»Top  view  Theta  values  for  quadrant  II 

Theta_tIII  =  []  ;  7»Top  view  Theta  values  for  quadrant  III 

Theta_tIV  =  []  ;  7»Top  view  Theta  values  for  quadrant  IV 

7«Side  view  quadrant  I  is  upper  left  corner  due  to 
7»  90  deg  theta  rotation 

Theta_sl  =  []  ;  7»Side  view  Theta  values  for  quadrant  I 

Theta_sII  =  []  ;  7»Side  view  Theta  values  for  quadrant  II 

Theta_sIII  =  []  ;  7»Side  view  Theta  values  for  quadrant  III 

Theta_sIV  =  []  ;  7»Side  view  Theta  values  for  quadrant  IV 

for  k  =  l:length(r) 

for  j  =  l:length(Xn) 

Theta_tII (j ,k)  =  atan2(MIDt_pos(j ,k) ,X(j)) ; 

7o  if  abs  (Theta_tII  ( j  ,k) )  ==  pi 

7o  Theta_tII  ( j  ,k)  =0; 

7o  end 

Theta_tIII ( j ,k)  =  atan2 (MIDt_neg( j ,k) ,X(j)) ; 

7o  if  abs  (Theta_tIII  ( j  ,k) )  ==  pi 

7o  Theta_tIII  ( j  ,k)  =0; 

7o  end 

end 

for  j  =  l:length(Xp) 

Theta_tl ( j ,k)  =  atan2(MIDt_pos(length(Xn)+j ,k) , . . . 
X(length(Xn)+j) ) ; 

Theta_tIV(j ,k)  =  atan2 (MIDt_neg (length (Xn) +j ,k) ,.. . 
X(length(Xn)+j) ) ; 

end 

for  j  =  l:length(Yn) 

Theta_sIV(j ,k)  =  atan2(Y(j) ,MIDs_pos(j ,k)) ; 

Theta_sIII ( j ,k)  =  atan2 (Y ( j ) ,MIDs_neg( j ,k) ) ; 

end 

for  j  =  l:length(Yp) 

Theta_sl (j ,k)  =  atan2(Y(length(Yn)+j) , . . . 
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MIDs_pos (length(Yn)+j ,k) ) ; 

Theta_sII ( j ,k)  =  atan2 (Y (length (Yn)+j ),..  . 
MIDs_neg(length(Yn)+j ,k) ) ; 

end 

end 


"/(Hemisphere  Theta  Matricies 
Theta_tN  =  [Theta_tII;  Theta_tl] ; 
Theta_tS  =  [Theta_tIII;  Theta_tIV] ; 
Theta_sE  =  [Theta_sIII;  Theta_sII] ; 
Theta_sW  =  [Theta_sIV;  Theta_sl] ; 


%Top  view  upper  half 
%Top  view  lower  half 
%Side  view  left  half 
%Side  view  right  half 


“/oSUBMATRIX  GENERATION 
/(Radially  dependent  submatrix 
Ml  =  [T;S]  ; 

[a,b]  =  size (Ml) ; 


/(Cosine  dependent  submatrix 

topC  =  [cos(Theta_tN)+cos(Theta_tS)]  ; 

sideC  =  [cos(Theta_sE)+cos(Theta_sW)] ; 

fullC  =  [topC; sideC] ; 

for  j  =  1 : a 

for  k  =  l:length(r) 

M2 ( j , k)  =  (Ml ( j ,k) /2) . *f ullC ( j ,k) ; 

end 

end 


°/0Sine  dependent  submatrix 

topS  =  [sin(Theta_tN)+sin(Theta_tS)]  ; 

sideS  =  [sin(Theta_sE)+sin(Theta_sW)] ; 

fullS  =  [topS ; sideS] ; 

for  j  =  1 : a 

for  k  =  l:length(r) 

M3 ( j , k)  =  (Ml (j ,k)/2) . *fullS(j ,k) ; 

end 

end 


%FULL  GEOMETRIC  MATRIX 
M  =  zeros (a,3*length(r) ) ; 


p  =  1; 
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for  j  =  l:length(r) 

M(: ,p)  =  MIC : ,j); 

M(: ,p+l)  =  M2 ( : , j ) ; 

M(: ,p+2)  =  M3 ( : , j ) ; 

P  =  P+3; 

end 

7.AXI SYMMETRIC  1ST  SHELL 

M  ( :  ,2:3)  =  []  ;  °/0force  axi  symmetric  first  shell 

[number_E,number_e] =size (M)  ; 

%axisymmetric ,  cosine,  and  sine  component  indicies 

10  =  2 : 3 :  number_e ;  °/0axi symmetric  indicies 

11  =  3 : 3 : number_e ;  °/0cosine  indicies 

12  =  4 : 3 :  number_e ;  °/0sine  indicies 

10=  [1,10];  °/„axi  symmetric  for  truncated  M  matrix 

°/„LEAST  SQUARES  FIT 
erraw  =  zeros (number_e , 1) ; 

Ecal  =  zeros (number_E, 1) ; 

erraw  =  M\signal; 

/(Constraint  on  Cosine  and  Sine  components 
eraxi  =  erraw (10) ; 

ervar  =  2*sqrt (erraw(Il) . ~2+erraw(I2) . ~2) ; 
for  j  =  l:length(I2) 

if  eraxi (j+1)  <  ervar(j) 
eraxi(j+l)  =  ervar(j); 

end 

er(I0(j+l))  =  eraxi(j+l); 

end 

er(l)  =  eraxi(l) ; 
er(Il)  =  erraw (II); 
er(I2)  =  erraw(I2); 
erT  =  erJ  ; 


